基于陀螺飞轮的航天器角速率实时线性化测量方法与流程

本发明是基于陀螺飞轮的航天器角速率实时线性化测量方法，具体涉及惯性导航领域。

背景技术：

陀螺飞轮可作为航天器的执行机构，为航天器提供三维姿态控制力矩，同时，与传统的动力调谐陀螺仪类似，陀螺飞轮可作为敏感器，实现航天器的两维角速率测量。为了实现陀螺飞轮的角速率测量功能，需提出切实可行的基于陀螺飞轮的航天器角速率测量方程。现有的利用陀螺飞轮实现两维航天器角速率测量方法，往往是利用广义坐标表示的陀螺飞轮系统方程，考虑陀螺飞轮电机转速恒定时，对该系统方程在零倾侧角位置处进行线性化，整理得到角速率测量方程，利用该思路得到的测量方程存在问题如下：

(1)现有方法根据动力调谐陀螺仪原理，对用广义坐标表示的陀螺飞轮动力学方程进行线性化时，线性化工作点取为零倾侧位置，而不同于动力调谐陀螺仪，陀螺飞轮为输出径向控制力矩，转子的工作倾侧角最大可达7°，此时，若采用现有测量方法进行角速率测量，同样将产生较大的测量误差；综上，现有基于陀螺飞轮的航天器角速率测量方法仅适用于陀螺飞轮电机转速恒定、且工作倾侧角在零倾侧位置小范围内的工作条件下，而这一工作条件严重限制了陀螺飞轮作为执行器输出三维控制力矩功能，从而无法充分发挥陀螺飞轮的结构与功能优势，大大制约了陀螺飞轮在微小航天器领域的应用；

(2)在前期研究中，根据完整陀螺飞轮非线性动力学方程，推导得到了一种非线性的陀螺飞轮系统测量方程，但是，由于测量方程本身含有大量的非线性项，导致测量方程的计算复杂度较高，从而对实际工程应用中软硬件实现实时性提出了较大的挑战。

技术实现要素：

针对现有技术的上述缺陷，本发明要解决现有基于陀螺飞轮的航天器角速率测量方法中测量精度受制于转子倾角较小及实时性差的问题，使陀螺飞轮力矩输出功能与角速率测量功能得到充分的应用。包括以下步骤：

步骤一、建立陀螺飞轮系统的运动学方程；

步骤二、建立陀螺飞轮系统的动力学方程；

步骤二一、构建陀螺飞轮系统的能量方程；

步骤二二、利用第二类拉格朗日法，建立陀螺飞轮系统的动力学方程；

步骤三、陀螺飞轮非线性动力学方程坐标变换

步骤三一、建立广义坐标与转子倾侧传感器可测倾角之间关系；

步骤三二、建立以壳体坐标表征的陀螺飞轮非线性动力学方程

步骤四、Lypapunov线性化陀螺飞轮动力学方程

步骤四一、Lypapunov线性化陀螺飞轮动力学方程

步骤四二、陀螺飞轮线性动力学约束方程建立

步骤五、基于陀螺飞轮的实时线性测量方程实现角速率测量

本发明有益效果：

本发明提出了一种基于陀螺飞轮的航天器角速率的实时线性化测量方法，该方法克服了零倾侧线性化所引起的不适于大倾侧工作状态的缺陷，可用于陀螺飞轮大倾角工作状态的航天器两维角速率测量；而且，本发明对非线性动力学方程进行了合理的线性化运算，避免了由于非线性而导致的测量方程计算复杂，实时性差的问题，利于工程实现。

附图说明

图1是本发明的方法流程图；

图2是本发明所用陀螺飞轮装置的简易结构图；

图3是本发明和现有方法所实现航天器x轴角速率ω_bx测量的对比仿真图；

图4是本发明和现有方法所实现航天器y轴角速率ω_by测量的对比仿真图。

具体实施方式

具体实施方式一：本实施方式基于陀螺飞轮的航天器角速率实时线性化测量方法按照以下步骤来实现：

步骤一、建立陀螺飞轮系统的运动学方程；

步骤二、建立陀螺飞轮系统的动力学方程；

步骤三、陀螺飞轮非线性动力学方程坐标变换；

步骤四、Lypapunov线性化陀螺飞轮动力学方程；

步骤五、基于陀螺飞轮的实时线性测量方程实现角速率测量。

具体实施方式二：本实施方式与具体实施方式一不同的是：其特征在于，所述的步骤一建立陀螺飞轮系统的运动学方程按照以下步骤实现：

陀螺飞轮系统核心构件由壳体(1)、电机轴(2)、平衡环(3)、转子(4)以及外扭杆(5)、内扭杆(6)组成，其中电机轴(2)通过一对内扭杆(6)与平衡环(3)内侧相连，平衡环(3)外侧通过一对外扭杆(5)与转子(4)内侧相连，内扭杆(6)与外扭杆(5)保持正交；

设陀螺飞轮系统的壳体、电机轴、平衡环、转子的体坐标系分别为Ox_cy_cz_c、Ox_my_mz_m、Ox_gy_gz_g、Ox_ry_rz_r，壳体、电机轴、平衡环、转子的体坐标系间的相对角位置关系如下：将Ox_cy_cz_c绕z_c轴旋转θ_z，得到Ox_my_mz_m，再将Ox_my_mz_m绕x_m旋转θ_x，得到Ox_gy_gz_g，最终将Ox_gy_gz_g绕y_g旋转θ_y，得到Ox_ry_rz_r；

其中，θ_z表示电机轴的转角，θ_x表示与平衡环连接的内扭杆转角，θ_y表示与平衡环连接的外扭杆转角；

据陀螺飞轮系统壳体与航天器固联，设定在陀螺飞轮系统壳体体坐标系下航天器相对惯性空间的角速度ω_b表示为ω_b＝[ω_bx ω_by ω_bz]^T：

则根据壳体、电机轴、平衡环、转子的体坐标系之间的相对位置关系，可分别求得电机轴、平衡环、转子在各自体坐标系下的转动角速度分别如下(1)(2)(3)所示：

其中，与分别表示转角θ_i的余弦值cosθ_i和正弦值sinθ_i；

i＝x，y，z；ω_m,ω_g,ω_r分别表示电机轴、平衡环、转子相对惯性空间的角速度。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：其特征在于，所述的步骤二建立陀螺飞轮系统的动力学方程按照以下步骤实现：

步骤二一、构建陀螺飞轮系统的能量方程；

陀螺飞轮系统的动能T由电机轴、平衡环及转子三部分转动动能组成，如公式(4)动能T的广义速度二次型表示为：

其中，I_mi，I_gi，I_ri分别表示为电机轴，平衡环及转子在惯性主轴上的转动惯量，i＝x，y，z；陀螺飞轮系统的势能V由陀螺飞轮系统内、外两扭杆的弹性形变引起，如公式(5)所示：

其中，K_x，K_y分别表示内扭杆、外扭杆的抗扭刚度；

陀螺飞轮系统的拉格朗日函数L即能量方程如公式(6)所示：

步骤二二、利用第二类拉格朗日法，建立陀螺飞轮系统的动力学方程；

陀螺飞轮系统的拉格朗日方程，如公式(7)所示：

其中，T_gx，T_gy，T_gz分别表示广义控制力矩分别在内扭杆、外扭杆及电机轴方向的投影；

C_x，C_y分别表示内、外两对扭杆的阻尼系数；T_fz为电机轴轴承产生的摩擦力矩；

由于公式(7)中第三式摩擦力矩T_fz的难以精确辨识，且公式(7)中第三式所表征的电机轴运动与公式(7)中第一式、第二式分别表征的转子沿两径向轴运动之间耦合微弱，仅需考虑公式(7)第一式、第二式进行动力学建模，将公式(6)拉格朗日函数带入至公式(7)第一式、第二式，得到陀螺飞轮系统动力学方程如公式(8)所示：

其中，

N₃＝I_ry cos θ_x sinθ_z N₄＝-I_ry cosθ_x cosθ_z

I_l＝I_gx+I_rx cos²θ_y+I_rz sin²θ_y I₂＝I_gz-I_gy-I_ry+I_rxsin²θ_y+I_rz cos²θ_y；

其中，设定动力学方程中航天器z轴方向的角速率ω_bz＝0，外部控制力矩T_gx，T_gy分别表示为：

其中，k_tx，k_ty均表示二维力矩器的标度因数；i_x，i_y均表示二维力矩器的电流。

其它步骤及参数与具体实施方式一至二之一相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：其特征在于，所述的步骤三陀螺飞轮非线性动力学方程坐标变换按照以下步骤实现：

步骤三一、建立广义坐标θ_x，θ_y与转子倾侧传感器可测倾角之间关系；

由于广义坐标θ_x，θ_y分别表示的为陀螺飞轮内扭杆和外扭杆的转角，而实际陀螺飞轮装置中的倾侧角传感器实现的为壳体系下的转子倾侧角的直接测量，无法直接测量得到θ_x，θ_y的测量值，而且本发明中的Lyapunov线性化也是针对壳体下转子的工作倾角进行的，因而，需定义壳体系Ox_cy_cz_c下的转子沿Ox_c和Oy_c的倾侧角φ_x，φ_y，根据电机轴、内扭杆、外扭杆以及壳体之间的相对位置关系，可得广义坐标θ_x，θ_y与转子倾侧传感器可测倾角φ_x，φ_y之间关系，如下：

相应地，根据方程(10)可得各阶导数之间的关系，如下：

步骤三二、壳体坐标表征的陀螺飞轮非线性动力学方程；

假设转子在横轴方向的转动惯量I_rx，I_ry相等且大小为I_ri，平衡环在横轴方向的转动惯量I_gx，I_gy相等且大小为I_gt，则满足如下关系：

将方程(10)-(13)带入至方程(8)中，并整理可得如下在壳体系下陀螺飞轮的非线性动力学方程：

其中，

系数矩阵元素具体如下：

I_re＝I_rs-I_ri

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：其特征在于，所述的步骤四Lypapunov线性化陀螺飞轮动力学方程按照以下步骤来实现：

步骤四一、Lypapunov线性化陀螺飞轮动力学方程

对陀螺飞轮系统，给定标称工作点x_d，ω_b，分别如下式所示：

其中，φ_xd，φ_yd分别为欲线性化的转子沿径向两轴的标称工作倾角，ω_bd为两轴外界航天器转动角速度矢量；

对壳体系下陀螺飞轮系统动力学方程(13)进行Lyapunov线性化，可得经李雅普诺夫线性化后的陀螺飞轮线性动力学方程如下：

其中，

I_ld＝I_gt+L_rt+I_re sin²θ_yd I_2d＝I_gs-I_gt+I_re cos²θ_yd I_re＝I_rs-I_rt

表示f_ci分别对φ_j，ω_bj，ω_bj进行求偏导；

T_cd＝[T_cxd T_cyd]^T为标称的两维力矩器控制力矩输入；

进一步计算整理方程(15)后可得Lyapunov线性化的陀螺飞轮动力学方程，如下(16)所示：

其中，

步骤四二、陀螺飞轮线性动力学约束方程建立

由于标称工作点(x_d，T_cd，ω_bd)同时满足方程(13)，因而，结合方程(14)，有如下方程(17)成立：

从而得到以标称工作倾角x_d表征的标称控制输入力矩T_cd，即陀螺飞轮线性动力学约束方程，如方程(18)所示：

由此，可得由方程(16)(18)构成的陀螺飞轮Lyapunov线性化动力学方程的完整表征形式。

其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：其特征在于，所述的步骤五基于陀螺飞轮实时线性测量方程实现角速率测量；

整理Lyapunov线性动力学方程(16)，可得陀螺飞轮的实时线性测量方程如下(19)所示：

其中，x_d＝[φ_xd φ_yd]^T为陀螺飞轮的两维输入倾侧角指令，即期望的陀螺飞轮标称倾侧工作状态，(φ_x，φ_y)为陀螺飞轮的实时工作倾侧角，可利用倾侧角传感器直接可测；同时，将倾侧角传感器测得到的(φ_x，φ_y)分别进行一次、二次差分可得到相对应的实时倾侧角速度和倾侧角加速度信息；(T_cx，T_cy)为力矩器输出力矩，如前一致，可用(k_txi_x，k_tyi_y)进行表征，其中两维力矩器电流(i_x，i_y)传感器可测，标度因子(k_tx，k_ty)可通过离线标定得到；T_cd＝[T_cxd T_cyd]^T为两维力矩器标称控制力矩，服从约束方程(18)，计算易得；测量方程(19)形式上为一阶微分方程，在实际工程中，对方程(19)进行离散化后，可采用递推方式实时求解两维航天器的角速率(ω_bx，ω_by)；

进一步，考虑搭载有陀螺飞轮实现角速率测量的航天器带宽很低，则实时线性测量方程方程(19)，可简化为：

该方程为标准的代数方程，结合传感器信息及系统参数，同样可实时实现基于陀螺飞轮的航天器角速率的实时线性测量；

其它步骤及参数与具体实施方式一至四之一相同。

实施例

为验证本发明所提方法的有效性，以一组陀螺飞轮系统参数为例，搭建仿真平台，对本发明所提方法以及现有的零倾侧线性化测量方法，进行对比仿真。陀螺飞轮的系统参数如下：

陀螺飞轮系统的转子在转子体坐标系下和平衡环在平衡环体坐标系下的沿惯性主轴方向的惯量参数分别为：

I_rt＝3.458×10^-3kg·m² I_gt＝1.2758×10^-5kg·m²

I_rs＝6.402×10^-3kg·m² I_gs＝1.8047×10^-5kg·m²

内、外扭杆的阻尼系数(c_x，c_y)均设为零，抗扭刚度分别选取为：

K_x＝0.092N·m/rad，K_y＝0.092N·m/rad；

电机转速为：

转子两维倾侧角(φ_xd，φ_yd)指令分别给定为：

φ_xd＝4·sin(2π·0.05·t)度

φ_xd＝0度

航天器的转动惯量J为：J＝diag[50 50 50]kg·m²，设定航天器初始角速度为零；

图3、图4分别实现的是利用陀螺飞轮实现的x轴和y轴航天器的角速率(ω_bx，ω_by)测量效果图；图3、图4中实线表示的是本发明提出的基于陀螺飞轮实时线性化角速率测量值；虚线表示的为真实的航天器角速率值；点线表示的是现有零倾侧线性测量方程的角速率测量值；从图3、图4的仿真结果可以看出，本发明提出的实时线性化角速率测量方法能够较为精确地实现航天器两维角速率的测量，而现有的零倾侧线性化测量方法在倾侧角为零附近时，可以实现准确的角速率测量，但是，如果陀螺飞轮转子工作倾角远大于零度，则该方法的测量精度急剧下降。本发明大大提高了陀螺飞轮在较大工作倾角下进行两维航天器角速率测量的能力，对陀螺飞轮同时作为执行器和敏感器在航天器姿态控制系统中的广泛应用具有重要意义。