CN103926840B - 一种主动抑制太阳帆板挠性振动的方法 - Google Patents

Abstract

一种主动抑制太阳帆板挠性振动的方法，该方法以正余弦细分驱动的步进电机为执行机构，建立卫星系统动力学模型，对系统动力学模型进行线性化化简后，计算系统的振动频率和阻尼比，设计输入成型器来对原始帆板转速指令进行整型，使得系统在完成帆板转速连续驱动的同时，还有效地抑制了帆板的挠性振动，从而降低了帆板挠性振动对载荷成像质量的影响，且对模型参数变化和设计参数不准具有较强的鲁棒性，相对于采用压电陶瓷材料和阻尼器等其它太阳帆板挠性振动控制方法，本发明不需要添加额外的硬件设备，工程实现比较简单。

Description

一种主动抑制太阳帆板挠性振动的方法

技术领域

本发明涉及一种主动抑制太阳帆板挠性振动的方法，属于挠性航天器控制技术领域，用于太阳帆板转速驱动切换过程中的挠性振动抑制控制，尤其是对低频和甚低频挠性振动的抑制。

背景技术

为了跟踪太阳以保证能源供应，目前卫星上普遍安装了太阳帆板及其驱动机构，但由于帆板的刚度不足，加上空间环境的弱阻尼特点，帆板的挠性振动一旦被激起后很难快速衰减。以往由于光学或SAR(合成孔径雷达)载荷的分辨率不高，导致对卫星姿态指向精度和稳定度等性能指标要求比较低，因此帆板挠性振动对载荷性能的影响并不明显。但是，随着载荷技术的发展，图像分辨率得到了很大的提高，相应地对姿态控制精度提出了越来越高的指标，稳定度要求达到了10-4°/s甚至更高的量级。在此背景下，太阳帆板挠性振动对载荷成像质量的影响就立马凸现出来，成为了制约遥感卫星进一步提升姿态控制性能的关键技术问题，因此，采用主动抑制方法来降低太阳帆板的挠性振动成了当前提高遥感卫星姿态控制精度的一个重要课题。

文献1“斯祝华，刘一武，黎康.太阳帆板驱动装置建模及其驱动控制研究[J].空间控制技术与应用，2010，36(2)：13～19”建立了基于正余弦细分驱动步进电机的太阳帆板驱动模型，并讨论了其驱动控制方式，但并未考虑帆板的挠性振动因素。

文献2“梁小光，丁竹生，焦映厚.卫星太阳翼阻尼器参数选定方法[J].哈尔滨工业大学学报，2011，43(7)：71～75”和文献3“胡庆雷，马广富.挠性航天器姿态机动控制的主动振动抑制[J].振动工程学报，18(3)：375～380”分别采用阻尼器和压电陶瓷材料对太阳帆板挠性进行主动抑制控制，这两种控制方式虽然能够降低帆板的挠性振动，但是需要增加敏感器和做动器等复杂硬件设备，星上工程实现尚未成熟。

除了文献2和3提出的挠性振动控制方法之外，其它挠性振动控制方法也得到了研究，其中，输入成型以其简单易用性、对模型参数摄动的鲁棒性等优良品质，逐渐成为目前主动振动抑制领域的主流方法。输入成型是指将脉冲序列与期望输入指令卷积形成的新指令作为系统输入的信号处理技术，其中脉冲序列称为输入成型器，脉冲序列的作用时间和幅值根据系统的模态频率和阻尼比得到，期望输入指令的确定仅考虑刚体系统运动的要求，整型后的指令不仅可以满足刚体的运动要求，还可以有效抑制系统中挠性附件的振动。

文献4“Michael J.Doherty，Robert J.Tolson.输入成型减少太阳帆板结构振动[C].美国航空航天局报告，No：19980232013”和文献5“原劲鹏，刘建功，杨雷.步进电机驱动柔性负载的一种振动抑制控制策略[J].空间控制技术与应用，2008，34(6)：34～38”虽然考虑了负载的挠性因素，采用输入成型技术实现了对挠性体的驱动定位和振动抑制，但是在它们的研究中，太阳帆板采用了斩波恒流驱动的步进电机，且应用背景仅限于挠性体的转角定位，转角到位之后电机停止驱动，这与目前工程上的太阳帆板采用正余弦细分驱动步进电机以及帆板需要连续驱动以跟踪太阳的实际工况尚存在差异。

发明内容

本发明提供一种主动抑制太阳帆板挠性振动的方法，可以有效抑制太阳帆板转速切换过程中的挠性振动，且对模型参数变化和设计参数不准具有比较强的鲁棒性，不需要添加额外的硬件设备，工程实现比较简单，具有较高的工程应用价值。

为了达到上述目的，本发明提供一种主动抑制太阳帆板挠性振动的方法，该方法包含以下步骤：

步骤1、建立太阳帆板挠性振动系统动力学模型并计算系统的振动模态参数；

步骤2、获得ZVD成型器；

步骤3、星载计算机利用ZVD综合成型器对太阳帆板的转速指令进行整形后，驱动步进电机。

所述的步骤1包含以下步骤：

步骤1.1、建立完整的太阳帆板挠性振动系统动力学模型：

所述的太阳帆板挠性振动系统包含卫星平台、一个太阳帆板和一个步进驱动电机，其振动系统动力学模型为：

$I\stackrel{\·}{\mathrm{\ω}}+{\mathrm{\ω}}^{\×}(\mathrm{I\ω}+h_w+R\stackrel{\·}{\mathrm{\Θ}}+F_s\stackrel{\·}{\mathrm{\η}})+R\stackrel{\·\·}{\mathrm{\Θ}}+F_b\stackrel{\·\·}{\mathrm{\η}}=T_d+T_c---(1-a)$

$J\stackrel{\·\·}{\mathrm{\Θ}}+R^T\stackrel{\·}{\mathrm{\ω}}+F_a\stackrel{\·\·}{\mathrm{\η}}+{\stackrel{\·}{\mathrm{\Θ}}}^{\×}(R^T\mathrm{\ω}+J\stackrel{\·}{\mathrm{\Θ}}+F_a\stackrel{\·}{\mathrm{\η}})=T_m---(1-b)$

$\stackrel{\·\·}{\mathrm{\η}}+2\mathrm{\ξ\Ω}\stackrel{\·}{\mathrm{\η}}+{\mathrm{\Ω}}^2\mathrm{\η}+F_b^T\stackrel{\·}{\mathrm{\ω}}+F_a^T\stackrel{\·\·}{\mathrm{\Θ}}=0---(1-c)$

其中，上标“^×”表示叉乘运算，ω∈R^3×1为卫星相对惯性系的角速度，R为实数集，为帆板转角，η∈R^n×1为帆板挠性振动模态坐标，h^w∈R^3×1为飞轮角动量，T_d∈R^3×1为卫星受到的环境力矩，T_c∈R^3×1为飞轮控制力矩，T_m∈R^3×1为步进电机驱动力矩，I∈R^3×3为卫星的惯量矩阵，J∈R^3×3为帆板的惯量矩阵，R∈R^3×3为帆板旋转与星体转动的耦合系数，F_b∈R^3×n为帆板振动运动与星体转动运动的耦合系数，F_a∈R^3×n为帆板振动运动与帆板驱动运动的耦合系数，ξ∈R^n×n为帆板挠性模态阻尼比，Ω∈R^n×n为帆板挠性模态频率；

步骤1.2、对系统动力学模型进行线性化化简：

令为3-1-2转序(即ZXY转序)下的卫星姿态，其中，是卫星的俯仰角，θ是卫星的滚转角，ψ是卫星的偏航角；

线性化化简的基础为卫星稳态飞行时的姿态角为小角度，因此卫星的惯性角速度和角加速度可分别化简为：

其中，ω₀为卫星的轨道角速度，A＝[0 0 -ω₀；0 0 0；ω₀ 0 0；]；

卫星姿态采用PD(比例-微分)控制，即：

$T_c=K_p({\mathrm{\Φ}}_0-\mathrm{\Φ})+K_d({\stackrel{\·}{\mathrm{\Φ}}}_0-\stackrel{\·}{\mathrm{\Φ}})---\left(4\right)$

其中，Φ₀为指令姿态，为指令姿态角速度， diag为对角阵；

太阳帆板为一维驱动，其驱动轴垂直于轨道面，指令转角 ${\mathrm{\Θ}}_r={\left[\begin{array}{ccc}0& {\mathrm{\θ}}_r& 0\end{array}\right]}^T,$ 其驱动力矩为：

$T_m=K{\mathrm{\Θ}}_r-B\stackrel{\·}{\mathrm{\Θ}}-\mathrm{K\Θ}---\left(5\right)$

其中，K＝[0 0 0；0 k 0；0 0 0；]，k为与驱动电机齿数、电流幅值等有关的系数，B＝[0 0 0；0 B 0；0 0 0；]，B为粘滞阻尼系数；

将公式(2)～(5)代入公式(1)并简化，得到系统的线性动力学模型为：

$\left\{\begin{array}{c}I\stackrel{\·\·}{\mathrm{\Φ}}+(\mathrm{IA}+C+K_d)\stackrel{\·}{\mathrm{\Φ}}+(D+K_p)\mathrm{\Φ}+R\stackrel{\·\·}{\mathrm{\Θ}}+F_b\stackrel{\·\·}{\mathrm{\η}}=T_{d1}+K_p{\mathrm{\Φ}}_0+K_d{\stackrel{\·}{\mathrm{\Φ}}}_0\\ J\stackrel{\·\·}{\mathrm{\Θ}}+R^T\stackrel{\·\·}{\mathrm{\Φ}}+R^{T}A\stackrel{\·}{\mathrm{\Φ}}+F_a\stackrel{\·\·}{\mathrm{\η}}+b\stackrel{\·}{\mathrm{\Θ}}+\mathrm{k\Θ}=K{\mathrm{\Θ}}_r+T_{d2}\\ \stackrel{\·\·}{\mathrm{\η}}+2\mathrm{\ξ\Ω}\stackrel{\·}{\mathrm{\η}}+{\mathrm{\Ω}}^2\mathrm{\η}+F_b^T\stackrel{\·\·}{\mathrm{\Φ}}+F_b^T\text{A}\stackrel{\·}{\mathrm{\Φ}}+F_a^T\stackrel{\·\·}{\mathrm{\Θ}}=0\end{array}\right.---\left(6\right)$

步骤1.3、计算系统的振动频率和阻尼比：

设状态变量 $X={\left[\begin{array}{cccccc}\mathrm{\Φ}& \stackrel{\·}{\mathrm{\Φ}}& \mathrm{\Θ}& \stackrel{\·}{\mathrm{\Θ}}& \mathrm{\η}& \stackrel{\·}{\mathrm{\η}}\end{array}\right]}^T,$ 则公式(6)可以写为如下状态方程：

$\stackrel{\·}{X}=A_1^{-1}{A}_{2}X+A_1^{-1}V---\left(7\right)$

系统矩阵具有2×(n+6)个复数特征值λ_i，其中，i＝1～2×(n+6)，则系统第i个振动模态的振动频率ω_i和阻尼比ζ_i与特征值λ_i存在以下关系：

${\mathrm{\λ}}_i=-{\mathrm{\ζ}}_i{\mathrm{\ω}}_i\±\sqrt{{{\mathrm{\ζ}}_i}^2-1}---\left(8\right)$

根据公式(8)可计算得到系统的振动模态的振动频率为阻尼比为ζ_i＝|a/ω_i|，其中，a和b分别为λ_i的实部和虚部；

系统有不止一个振动频率及其对应的阻尼比，有n个振动模态就有n个振动模态，振动频率ω_i和阻尼比ζ_i只是其中的一个。

所述的步进驱动电机采用正余弦细分驱动的步进电机。

所述的步骤2包含以下步骤：

ZVD成型器是一个由三个子脉冲组成的脉冲序列，其数学表达式为：

$I\left(t\right)=\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{A}_i\mathrm{\δ}(t-t_i)---\left(9\right)$

其中，A_i和t_i分别为第i个脉冲的幅值和施加时刻；

定义如下二阶系统：

$G\left(s\right)=\frac{{\mathrm{\ω}}^2}{s^2+2\mathrm{\ζ\ωs}+{\mathrm{\ω}}^2}---\left(10\right)$

其中，ω和ζ分别为系统的振动频率(也即系统的自然频率)和系统的阻尼比，将系统对脉冲序列公式(9)的响应与对未成型脉冲信号的响应的幅值之比定义为振动比：

$V=e^{-\mathrm{\ζ\ω}{t}_N}\sqrt{{\left[\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_i{e}^{-\mathrm{\ζ\ω}{t}_i}\cos \left(\mathrm{\ω}\sqrt{1-{\mathrm{\ζ}}^2}{t}_i\right)\right]}^2+{[\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_i{e}^{-\mathrm{\ζ\ω}{t}_i}\sin \left(\mathrm{\ω}\sqrt{1-{\mathrm{\ζ}}^2}{t}_i\right)}^2}---\left(11\right)$

ZVD成型器要求当最后一个脉冲A_Nδ(t-t_N)作用完毕之后，系统的振动为零，即要求下式成立：

$\left\{\begin{array}{c}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_i{e}^{-\mathrm{\ζ\ω}{t}_i}\cos \left(\mathrm{\ω}\sqrt{1-{\mathrm{\ζ}}^2}{t}_i\right)=0\\ \underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_i{e}^{-\mathrm{\ζ\ω}{t}_i}\sin \left(\mathrm{\ω}\sqrt{1-{\mathrm{\ζ}}^2}{t}_i\right)=0\end{array}\right.---\left(12\right)$

由于公式(12)存在多值解，对此需要增加以下约束条件：

1、为了保证输入成型器是物理可实现的，要求t_i≥0，且第一个脉冲出现在零时刻，要求t₁＝0；

2、所有脉冲幅值之和必须等于原始信号幅值，要求

3、脉冲幅值可正可负，在此约束为A_i>0；

4、为提高成型器的鲁棒性，增加频率导数约束dV/dω＝0；

根据以上约束条件对公式(12)求解得到ZVD成型器为：

$\left[\begin{array}{c}A\\ t\end{array}\right]=\left[\begin{array}{ccc}{A}_1& A_2& A_3\\ t_1& t_2& t_3\end{array}\right]\left[\begin{array}{ccc}\frac{1}{1+2K+K^2}& \frac{2K}{1+2K+K^2}& \frac{K^2}{1+2K+K^2}\\ 0& \frac{T_d}{2}& T_d\end{array}\right]---\left(13\right)$

其中， $K=e^{-(\mathrm{\ζ\π}/\sqrt{1-{\mathrm{\ζ}}^2})},T_d=\frac{2\mathrm{\π}}{\mathrm{\ω}\sqrt{1-{\mathrm{\ζ}}^2}}.$

所述的步骤3包含以下步骤：

步骤3.1、星载计算机将所有振动模态的ZVD成型器卷积得到一个ZVD综合成型器；

步骤3.2、星载计算机将太阳帆板的原始转速指令送入ZVD综合成型器整型后，将获得的整形指令作为实际转速指令来驱动步进电机。

所述的主动抑制太阳帆板挠性振动的方法还包含步骤4，星载计算机在轨修正ZVD成型器。

所述的步骤4包含以下步骤：

步骤4.1、通过对星载计算机下传的卫星姿态数据进行辨识，得到包括系统振动频率和阻尼比在内的系统振动模态参数的辨识值；

步骤4.2、判断振动模态参数的理论值和真实值之间的差别，若振动模态参数的理论值和辨识值之间相差大于10％，则根据辨识结果重新计算ZVD成型器的参数，并通过在轨编程方式修正存储在星载计算机中的ZVD成型器参数；如果理论值和辨识值相差小于等于10％，则不进行星上修正。

本发明针对当前工程上太阳帆板驱动的实际情况，以太阳帆板连续驱动控制为应用背景，基于输入成型技术，提出了一种以正余弦电流细分型步进电机驱动的太阳帆板挠性振动抑制控制方法。该方法以正余弦细分驱动的步进电机为执行机构，根据卫星系统动力学模型的振动模态参数，设计输入成型器来对原始帆板转速指令进行整型，使得系统在完成帆板转速连续驱动的同时，还有效地抑制了帆板的挠性振动，从而降低了帆板挠性振动对载荷成像质量的影响，且对模型参数变化和设计参数不准具有较强的鲁棒性；此外，相对于采用压电陶瓷材料和阻尼器等其它太阳帆板挠性振动控制方法，本发明不需要添加额外的硬件设备，工程实现比较简单。综合而言，本发明具有明显的优点和良好的工程应用价值。

附图说明

图1为本发明的流程图。

图2为ZVD输入成型器在太阳帆板驱动上的实现示意图。

图3为帆板挠性振动抑制效果图。

具体实施方式

以下根据图1～图3，具体说明本发明的较佳实施例。

如图1所示，本发明提供一种主动抑制太阳帆板挠性振动的方法，本方法包含以下步骤：

步骤1、建立太阳帆板挠性振动系统动力学模型并计算系统的振动模态参数；

步骤1.1、建立完整的太阳帆板挠性振动系统动力学模型：

所述的太阳帆板挠性振动系统包含卫星平台、一个太阳帆板和一个步进驱动电机(所述的步进驱动电机采用正余弦细分驱动的步进电机)，其振动系统动力学模型为：

$I\stackrel{\·}{\mathrm{\ω}}+{\mathrm{\ω}}^{\×}(\mathrm{I\ω}+h_w+R\stackrel{\·}{\mathrm{\Θ}}+F_s\stackrel{\·}{\mathrm{\η}})+R\stackrel{\·\·}{\mathrm{\Θ}}+F_b\stackrel{\·\·}{\mathrm{\η}}=T_d+T_c---(1-a)$

$J\stackrel{\·\·}{\mathrm{\Θ}}+R^T\stackrel{\·}{\mathrm{\ω}}+F_a\stackrel{\·\·}{\mathrm{\η}}+{\stackrel{\·}{\mathrm{\Θ}}}^{\×}(R^T\mathrm{\ω}+J\stackrel{\·}{\mathrm{\Θ}}+F_a\stackrel{\·}{\mathrm{\η}})=T_m---(1-b)$

$\stackrel{\·\·}{\mathrm{\η}}+2\mathrm{\ξ\Ω}\stackrel{\·}{\mathrm{\η}}+{\mathrm{\Ω}}^2\mathrm{\η}+F_b^T\stackrel{\·}{\mathrm{\ω}}+F_a^T\stackrel{\·\·}{\mathrm{\Θ}}=0---(1-c)$

其中，上标“^×”表示叉乘运算，ω∈R^3×1为卫星相对惯性系的角速度，R为实数集，为帆板转角，η∈R^n×1为帆板挠性振动模态坐标，h_w∈R^3×1为飞轮角动量，T_d∈R^3×1为卫星受到的环境力矩，T_c∈R^3×1为飞轮控制力矩，T_m∈R^3×1为步进电机驱动力矩，I∈R^3×3为卫星的惯量矩阵，J∈R^3×3为帆板的惯量矩阵，R∈R^3×3为帆板旋转与星体转动的耦合系数，F_b∈R^3×n为帆板振动运动与星体转动运动的耦合系数，F_a∈R^3×n为帆板振动运动与帆板驱动运动的耦合系数，ξ∈R^n×n为帆板挠性模态阻尼比，Ω∈R^n×n为帆板挠性模态频率；

步骤1.2、对系统动力学模型进行线性化化简：

公式(1)给出的系统动力学模型是一个非线性方程，不便于分析系统特性和设计输入成型器，因此需要进行线性化化简；

令为3-1-2转序(即ZXY转序)下的卫星姿态，其中，是卫星的俯仰角，θ是卫星的滚转角，ψ是卫星的偏航角；

线性化化简的基础为卫星稳态飞行时的姿态角为小角度，因此卫星的惯性角速度和角加速度可分别化简为

其中，ω₀为卫星的轨道角速度，A＝[0 0 -ω₀；0 0 0；ω₀ 0 0；]；

卫星姿态采用PD(比例-微分)控制，即：

$T_c=K_p({\mathrm{\Φ}}_0-\mathrm{\Φ})+K_d({\stackrel{\·}{\mathrm{\Φ}}}_0-\stackrel{\·}{\mathrm{\Φ}})---\left(4\right)$

其中，Φ₀为指令姿态，为指令姿态角速度， diag为对角阵；

太阳帆板为一维驱动，其驱动轴垂直于轨道面，指令转角 ${\mathrm{\Θ}}_r={\left[\begin{array}{ccc}0& {\mathrm{\θ}}_r& 0\end{array}\right]}^T,$ 其驱动力矩为：

$T_m=K{\mathrm{\Θ}}_r-B\stackrel{\·}{\mathrm{\Θ}}-\mathrm{K\Θ}---\left(5\right)$

其中，K＝[0 0 0；0 k 0；0 0 0；]，k为与驱动电机齿数、电流幅值等有关的系数，B＝[0 0 0；0 B 0；0 0 0；]，B为粘滞阻尼系数；

将公式(2)～(5)代入公式(1)并简化，得到系统的线性动力学模型为：

$\left\{\begin{array}{c}I\stackrel{\·\·}{\mathrm{\Φ}}+(\mathrm{IA}+C+K_d)\stackrel{\·}{\mathrm{\Φ}}+(D+K_p)\mathrm{\Φ}+R\stackrel{\·\·}{\mathrm{\Θ}}+F_b\stackrel{\·\·}{\mathrm{\η}}=T_{d1}+K_p{\mathrm{\Φ}}_0+K_d{\stackrel{\·}{\mathrm{\Φ}}}_0\\ J\stackrel{\·\·}{\mathrm{\Θ}}+R^T\stackrel{\·\·}{\mathrm{\Φ}}+R^{T}A\stackrel{\·}{\mathrm{\Φ}}+F_a\stackrel{\·\·}{\mathrm{\η}}+b\stackrel{\·}{\mathrm{\Θ}}+\mathrm{k\Θ}=K{\mathrm{\Θ}}_r+T_{d2}\\ \stackrel{\·\·}{\mathrm{\η}}+2\mathrm{\ξ\Ω}\stackrel{\·}{\mathrm{\η}}+{\mathrm{\Ω}}^2\mathrm{\η}+F_b^T\stackrel{\·\·}{\mathrm{\Φ}}+F_b^T\text{A}\stackrel{\·}{\mathrm{\Φ}}+F_a^T\stackrel{\·\·}{\mathrm{\Θ}}=0\end{array}\right.---\left(6\right)$

步骤1.3、计算系统的振动模态的参数(振动频率和阻尼比)：

设状态变量 $X={\left[\begin{array}{cccccc}\mathrm{\Φ}& \stackrel{\·}{\mathrm{\Φ}}& \mathrm{\Θ}& \stackrel{\·}{\mathrm{\Θ}}& \mathrm{\η}& \stackrel{\·}{\mathrm{\η}}\end{array}\right]}^T,$ 则公式(6)可以写为如下状态方程：

$\stackrel{\·}{X}=A_1^{-1}{A}_{2}X+A_1^{-1}V---\left(7\right)$

系统矩阵具有2×(n+6)个复数特征值λ_i，其中，i＝1～2×(n+6)，则系统第i个振动模态的振动频率ω_i和阻尼比ζ_i与特征值λ_i存在以下关系：

${\mathrm{\λ}}_i=-{\mathrm{\ζ}}_i{\mathrm{\ω}}_i\±\sqrt{{{\mathrm{\ζ}}_i}^2-1}---\left(8\right)$

根据公式(8)可计算得到系统的振动模态的振动频率为阻尼比为ζ_i＝|a/ω_i|，其中，a和b分别为λ_i的实部和虚部；

系统有不止一个振动频率及其对应的阻尼比，有n个振动模态就有n个振动模态，振动频率ω_i和阻尼比ζ_i只是其中的一个；

步骤2、获得ZVD成型器；

本发明采用ZVD成型器，这种成型器的优点是对系统振动频率和阻尼比变化和估计不准具有比较强的鲁棒性，并且对原始指令的响应滞后比较短；

ZVD成型器是一个由三个子脉冲组成的脉冲序列，其数学表达式为：

$I\left(t\right)=\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{A}_i\mathrm{\δ}(t-t_i)---\left(9\right)$

其中，A_i和t_i分别为第i个脉冲的幅值和施加时刻；

定义如下二阶系统：

$G\left(s\right)=\frac{{\mathrm{\ω}}^2}{s^2+2\mathrm{\ζ\ωs}+{\mathrm{\ω}}^2}---\left(10\right)$

其中，ω和ζ分别为系统的振动频率(也即系统的自然频率)和系统的阻尼比，将系统对脉冲序列公式(9)的响应与对未成型脉冲信号的响应的幅值之比定义为振动比：

$V=e^{-\mathrm{\ζ\ω}{t}_N}\sqrt{{\left[\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_i{e}^{-\mathrm{\ζ\ω}{t}_i}\cos \left(\mathrm{\ω}\sqrt{1-{\mathrm{\ζ}}^2}{t}_i\right)\right]}^2+{[\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_i{e}^{-\mathrm{\ζ\ω}{t}_i}\sin \left(\mathrm{\ω}\sqrt{1-{\mathrm{\ζ}}^2}{t}_i\right)}^2}---\left(11\right)$

ZVD成型器要求当最后一个脉冲A_Nδ(t-t_N)作用完毕之后，系统的振动为零，即要求下式成立：

$\left\{\begin{array}{c}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_i{e}^{-\mathrm{\ζ\ω}{t}_i}\cos \left(\mathrm{\ω}\sqrt{1-{\mathrm{\ζ}}^2}{t}_i\right)=0\\ \underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_i{e}^{-\mathrm{\ζ\ω}{t}_i}\sin \left(\mathrm{\ω}\sqrt{1-{\mathrm{\ζ}}^2}{t}_i\right)=0\end{array}\right.---\left(12\right)$

由于公式(12)存在多值解，对此需要增加以下约束条件：

1、为了保证输入成型器是物理可实现的，要求t_i≥0，且第一个脉冲出现在零时刻，要求t₁＝0；

2、所有脉冲幅值之和必须等于原始信号幅值，要求

3、脉冲幅值可正可负，在此约束为A_i>0；

4、为提高成型器的鲁棒性，增加频率导数约束dV/dω＝0；

根据以上约束条件对公式(12)求解得到ZVD成型器为：

$\left[\begin{array}{c}A\\ t\end{array}\right]=\left[\begin{array}{ccc}{A}_1& A_2& A_3\\ t_1& t_2& t_3\end{array}\right]\left[\begin{array}{ccc}\frac{1}{1+2K+K^2}& \frac{2K}{1+2K+K^2}& \frac{K^2}{1+2K+K^2}\\ 0& \frac{T_d}{2}& T_d\end{array}\right]---\left(13\right)$

其中， $K=e^{-(\mathrm{\ζ\π}/\sqrt{1-{\mathrm{\ζ}}^2})},T_d=\frac{2\mathrm{\π}}{\mathrm{\ω}\sqrt{1-{\mathrm{\ζ}}^2}}.$

以下通过一个实施例1来说明ZVD成型器的建立：

设拟抑制太阳帆板前两阶振动模态，频率分别为0.05Hz和0.30Hz，阻尼比分别为0.005和0.006，把这些参数代入公式(13)，得到ZVD成型器分别为：

ZVD成型器1(抑制0.05Hz振动)：

$\left[\begin{array}{c}{A}_1\\ t_1/s\end{array}\right]=\left[\begin{array}{ccc}{A}_{11}& A_{12}& A_{13}\\ 0& t_{12}& t_{13}\end{array}\right]=\left[\begin{array}{ccc}0.2539& 0.5000& 0.2461\\ 0& 10.000& 20.000\end{array}\right]$

ZVD成型器2(抑制0.30Hz振动)：

$\left[\begin{array}{c}{A}_2\\ t_2/s\end{array}\right]=\left[\begin{array}{ccc}{A}_{21}& A_{22}& A_{23}\\ 0& t_{22}& t_{23}\end{array}\right]=\left[\begin{array}{ccc}0.2551& 0.4999& 0.2450\\ 0& 1.6934& 3.3868\end{array}\right]$

步骤3、星载计算机利用ZVD成型器对太阳帆板的转速指令进行整形后，驱动步进电机；

步骤3.1、星载计算机将所有振动模态的ZVD成型器卷积得到一个综合ZVD成型器；

针对拟抑制的挠性振动模态(模态参数根据步骤1理论计算得到)，分别单独设计各自的ZVD成型器(根据步骤2给出的方法建立)，然后把所有成型器卷积得到一个综合ZVD成型器；

仍以实施例1为例，令符号“*”表示卷积运算，则能够抑制帆板前两阶振动模态的ZVD综合成型器＝ZVD成型器1*ZVD成型器2，即：

$\begin{array}{c}\left[\begin{array}{c}A\\ t/s\end{array}\right]=\left[\begin{array}{ccc}{A}_{11}& A_{12}& A_{13}\\ 0& t_{12}& t_{13}\end{array}\right]*\left[\begin{array}{ccc}{A}_{21}& A_{22}& A_{23}\\ 0& t_{22}& t_{23}\end{array}\right]\\ =\left[\begin{array}{ccccccccc}{A}_{11}{A}_{21}& A_{11}{A}_{22}& A_{11}{A}_{23}& A_{12}{A}_{21}& A_{12}{A}_{22}& A_{12}{A}_{23}& A_{13}{A}_{21}& A_{13}{A}_{22}& A_{13}{A}_{23}\\ 0& t_{22}& t_{23}& t_{12}& t_{12}+t_{22}& t_{12}+t_{23}& t_{13}& t_{13}+t_{22}& t_{13}+t_{23}\end{array}\right]\\ =\left[\begin{array}{ccccccccc}0.065& 0.127& 0.062& 0.128& 0.250& 0.123& 0.063& 0.123& 0.060\\ 0& 1.693& 3.387& 10.000& 11.693& 13.387& 20.000& 21.693& 23.387\end{array}\right]\end{array}$

步骤3.2、星载计算机将太阳帆板的原始转速指令送入ZVD综合成型器整型后，将获得的整形指令作为实际转速指令来驱动步进电机；

如图2所示，为ZVD输入成型器在太阳帆板驱动上的实现示意图，图中，变量A_ij表示第i个ZVD成型器(抑制系统第i个模态振动)的第j个脉冲的幅值，变量t_ij表示第i个ZVD成型器的第j个脉冲的作用时间，单个ZVD成型器在帆板转速驱动中的实现是将原始转速指令按照公式(13)中的矩阵第一列参数所确定的比例分成3部分，然后在由公式(13)中的矩阵第二列参数所确定的时刻(t₁₁＝0，t₁₂，t₁₃)依次把各部分和前面已输出的部分相加，相加之和作为新的转速指令再输给下一个ZVD成型器或步进电机；

将ZVD综合成型器添加到步进电机的驱动系统中，得到的2个太阳帆板挠性振动模态坐标如图3所示，从图中可以看出，整型后的转速指令有效抑制了帆板的挠性振动，振动幅值降低了至少一个数量级；

步骤4、星载计算机在轨修正ZVD成型器；

步骤1通过系统动力学模型计算得到了系统的振动模态参数，由于各种因素，理论计算结果与真实结果难免存在一定偏差，导致ZVD成型器存在参数误差，针对这个问题，在卫星发射上天后可通过在轨编程的方式修正ZVD成型器的参数误差，具体包含如下步骤：

步骤4.1、通过对星载计算机下传的卫星姿态数据进行辨识，得到包括系统振动频率和阻尼比在内的系统振动模态参数的辨识值(辨识值可看做是真实值)；

步骤4.2、判断振动模态参数的理论值和真实值之间的差别，若振动模态参数的理论值和辨识值之间相差大于10％(根据成型器鲁棒性决定)，则根据辨识结果重新计算ZVD成型器的参数，并通过在轨编程方式修正存储在星载计算机中的ZVD成型器参数；如果理论值和辨识值相差在允许范围内，则不进行星上修正。

本发明针对当前工程上太阳帆板驱动的实际情况，以太阳帆板连续驱动控制为应用背景，基于输入成型技术，提出了一种以正余弦电流细分型步进电机驱动的太阳帆板挠性振动抑制控制方法。该方法以正余弦细分驱动的步进电机为执行机构，根据卫星系统动力学模型的振动模态参数，设计输入成型器来对原始帆板转速指令进行整型，使得系统在完成帆板转速连续驱动的同时，还有效地抑制了帆板的挠性振动，从而降低了帆板挠性振动对载荷成像质量的影响，且对模型参数变化和设计参数不准具有较强的鲁棒性；此外，相对于采用压电陶瓷材料和阻尼器等其它太阳帆板挠性振动控制方法，本发明不需要添加额外的硬件设备，工程实现比较简单。综合而言，本发明具有明显的优点和良好的工程应用价值。

尽管本发明的内容已经通过上述优选实施例作了详细介绍，但应当认识到上述的描述不应被认为是对本发明的限制。在本领域技术人员阅读了上述内容后，对于本发明的多种修改和替代都将是显而易见的。因此，本发明的保护范围应由所附的权利要求来限定。

Claims (2)

1.一种主动抑制太阳帆板挠性振动的方法，其特征在于，该方法包含以下步骤：

步骤1、建立太阳帆板挠性振动系统动力学模型并计算系统的振动模态参数；

步骤2、获得ZVD成型器；

步骤3、星载计算机利用ZVD综合成型器对太阳帆板的转速指令进行整形后，驱动步进电机；

所述的步骤1包含以下步骤：

步骤1.1、建立完整的太阳帆板挠性振动系统动力学模型：

所述的太阳帆板挠性振动系统包含卫星平台、一个太阳帆板和一个步进驱动电机，其振动系统动力学模型为：

$I\stackrel{\·}{\mathrm{\ω}}+{\mathrm{\ω}}^{\×}(I\mathrm{\ω}+h_w+R\stackrel{\·}{\mathrm{\Θ}}+F_s\stackrel{\·}{\mathrm{\η}})+R\stackrel{\·\·}{\mathrm{\Θ}}+F_b\stackrel{\·\·}{\mathrm{\η}}=T_d+T_c---(1-a)$

$J\stackrel{\·\·}{\mathrm{\Θ}}+R^T\stackrel{\·}{\mathrm{\ω}}+F_a\stackrel{\·\·}{\mathrm{\η}}+{\stackrel{\·}{\mathrm{\Θ}}}^{\×}(R^T\mathrm{\ω}+J\stackrel{\·}{\mathrm{\Θ}}+F_a\stackrel{\·}{\mathrm{\η}})=T_m---(1-b)$

$\stackrel{\·\·}{\mathrm{\η}}+2\mathrm{\ξ}\mathrm{\Ω}\stackrel{\·}{\mathrm{\η}}+{\mathrm{\Ω}}^2\mathrm{\η}+F_b^T\stackrel{\·}{\mathrm{\ω}}+F_a^T\stackrel{\·\·}{\mathrm{\Θ}}=0---(1-c)$

其中，上标“×”表示叉乘运算，ω∈R′^3×1为卫星相对惯性系的角速度，R′为实数集，Θ∈R′^3×1为帆板转角，η∈R′^n×1为帆板挠性振动模态坐标，h_w∈R′^3×1为飞轮角动量，T_d∈R′^3×1为卫星受到的环境力矩，T_c∈R′^3×1为飞轮控制力矩，T_m∈R′^3×1为步进电机驱动力矩，I∈R′^3×3为卫星的惯量矩阵，J∈R′^3×3为帆板的惯量矩阵，R∈R′^3×3为帆板旋转与星体转动的耦合系数，F_b∈R′^3×n为帆板振动运动与星体转动运动的耦合系数，F_a∈R′^3×n为帆板振动运动与帆板驱动运动的耦合系数，ξ∈R′^n×n为帆板挠性模态阻尼比，Ω∈R′^n×n为帆板挠性模态频率；

步骤1.2、对系统动力学模型进行线性化化简：

令为3-1-2转序，即ZXY转序，下的卫星姿态，其中，是卫星的俯仰角，θ是卫星的滚转角，ψ是卫星的偏航角；

线性化化简的基础为卫星稳态飞行时的姿态角为小角度，因此卫星的惯性角速度和角加速度分别化简为：

其中，ω₀为卫星的轨道角速度，A＝[0 0 -ω₀；0 0 0；ω₀ 0 0；]；

卫星姿态采用PD，即比例-微分控制，即：

$T_c=K_p({\mathrm{\Φ}}_0-\mathrm{\Φ})+K_d({\stackrel{\·}{\mathrm{\Φ}}}_0-\stackrel{\·}{\mathrm{\Φ}})---\left(4\right)$

其中，Φ₀为指令姿态，为指令姿态角速度，diag为对角阵；

太阳帆板为一维驱动，其驱动轴垂直于轨道面，指令转角Θ_r＝[0 θ_r 0]^T，其驱动力矩为：

$T_m={\mathrm{K\Θ}}_r-B\stackrel{\·}{\mathrm{\Θ}}-K\mathrm{\Θ}---\left(5\right)$

其中，K＝[0 0 0；0 k 0；0 0 0；]，k为与驱动电机齿数、电流幅值有关的系数，B＝[0 00；0 B 0；0 0 0；]，B为粘滞阻尼系数；

将公式(2)～(5)代入公式(1)并简化，得到系统的线性动力学模型为：

$\left\{\begin{array}{c}I\stackrel{\·\·}{\mathrm{\Φ}}+(IA+C+K_d)\stackrel{\·}{\mathrm{\Φ}}+(D+K_p)\mathrm{\Φ}+R\stackrel{\·\·}{\mathrm{\Θ}}+F_b\stackrel{\·\·}{\mathrm{\η}}=T_{d1}+K_p{\mathrm{\Φ}}_0+K_d{\stackrel{\·}{\mathrm{\Φ}}}_0\\ J\stackrel{\·\·}{\mathrm{\Θ}}+R^T\stackrel{\·\·}{\mathrm{\Φ}}+R^{T}A\stackrel{\·}{\mathrm{\Φ}}+F_a\stackrel{\·\·}{\mathrm{\η}}+B\stackrel{\·}{\mathrm{\Θ}}+K\mathrm{\Θ}={\mathrm{K\Θ}}_r+T_{d2}\\ \stackrel{\·\·}{\mathrm{\η}}+2\mathrm{\ξ}\mathrm{\Ω}\stackrel{\·}{\mathrm{\η}}+{\mathrm{\Ω}}^2\mathrm{\η}+F_b^T\stackrel{\·\·}{\mathrm{\Φ}}+F_b^{T}A\stackrel{\·}{\mathrm{\Φ}}+F_a^T\stackrel{\·\·}{\mathrm{\Θ}}=0\end{array}\right.---\left(6\right)$

步骤1.3、计算系统的振动频率和阻尼比：

设状态变量则公式(6)写为如下状态方程：

$\stackrel{\·}{X}=A_1^{-1}{A}_{2}X+A_1^{-1}V---\left(7\right)$

系统矩阵具有2×(n+6)个复数特征值λ_i，其中，i＝1～2×(n+6)，则系统第i个振动模态的振动频率ω_i和阻尼比ζ_i与特征值λ_i存在以下关系：

${\mathrm{\λ}}_i=-{\mathrm{\ζ}}_i{\mathrm{\ω}}_i\±\sqrt{{\mathrm{\ζ}}_i^2-1}---\left(8\right)$

根据公式(8)计算得到系统的振动模态的振动频率为阻尼比为ζ_i＝|a/ω_i|，其中，a和b分别为λ_i的实部和虚部；

系统有不止一个振动频率及其对应的阻尼比，有n个振动模态就有n个振动频率，振动频率ω_i和阻尼比ζ_i只是其中的一个；

所述的步骤2包含以下步骤：

ZVD成型器是一个由三个子脉冲组成的脉冲序列，其数学表达式为：

$I\left(t\right)=\underset{i=1}{\overset{3}{\Σ}}{A}_i\mathrm{\δ}(t-t_i)---\left(9\right)$

其中，A_i和t_i分别为第i个脉冲的幅值和施加时刻；

定义如下二阶系统：

$G\left(s\right)=\frac{{\mathrm{\ω}}^2}{s^2+2\mathrm{\ζ}\mathrm{\ω}s+{\mathrm{\ω}}^2}---\left(10\right)$

其中，ω和ζ分别为系统的振动频率，也即系统的自然频率和系统的阻尼比，将系统对脉冲序列公式(9)的响应与对未成型脉冲信号的响应的幅值之比定义为振动比：

$V=e^{-{\mathrm{\ζ\ωt}}_N}\sqrt{{\[\underset{i=1}{\overset{N}{\Σ}}{A}_i{e}^{-{\mathrm{\ζ\ωt}}_i}cos\left(\mathrm{\ω}\sqrt{1-{\mathrm{\ζ}}^2}{t}_i\right)\]}^2+{\[\underset{i=1}{\overset{N}{\Σ}}{A}_i{e}^{-{\mathrm{\ζ\ωt}}_i}sin\left(\mathrm{\ω}\sqrt{1-{\mathrm{\ζ}}^2}{t}_i\right)\]}^2}---\left(11\right)$

ZVD成型器要求当最后一个脉冲A_Nδ(t-t_N)作用完毕之后，系统的振动为零，即要求下式成立：

$\left\{\begin{array}{c}\underset{i=1}{\overset{N}{\Σ}}{A}_i{e}^{-{\mathrm{\ζ\ωt}}_i}cos\left(\mathrm{\ω}\sqrt{1-{\mathrm{\ζ}}^2}{t}_i\right)=0\\ \underset{i=1}{\overset{N}{\Σ}}{A}_i{e}^{-{\mathrm{\ζ\ωt}}_i}sin\left(\mathrm{\ω}\sqrt{1-{\mathrm{\ζ}}^2}{t}_i\right)=0\end{array}\right.---\left(12\right)$

由于公式(12)存在多值解，对此需要增加以下约束条件：

a、为了保证输入成型器是物理可实现的，要求t_i≥0，且第一个脉冲出现在零时刻，要求t₁＝0；

b、所有脉冲幅值之和必须等于原始信号幅值，要求

c、脉冲幅值约束为A_i>0；

d、为提高成型器的鲁棒性，增加频率导数约束dV/dω＝0；

根据以上约束条件对公式(12)求解得到ZVD成型器为：

$\left[\begin{array}{c}A\\ t\end{array}\right]=\left[\begin{array}{ccc}{A}_1& A_2& A_3\\ t_1& t_2& t_3\end{array}\right]=\left[\begin{array}{ccc}\frac{1}{1+2K+K^2}& \frac{2K}{1+2K+K^2}& \frac{K^2}{1+2K+K^2}\\ 0& \frac{T_d}{2}& T_d\end{array}\right]---\left(13\right)$

其中，

所述的步骤3包含以下步骤：

步骤3.1、星载计算机将所有振动模态的ZVD成型器卷积得到一个ZVD综合成型器；

步骤3.2、星载计算机将太阳帆板的原始转速指令送入ZVD综合成型器整型后，将获得的整形指令作为实际转速指令来驱动步进电机；

所述的主动抑制太阳帆板挠性振动的方法还包含步骤4，星载计算机在轨修正ZVD成型器；

所述的步骤4包含以下步骤：

步骤4.1、通过对星载计算机下传的卫星姿态数据进行辨识，得到包括系统振动频率和阻尼比在内的系统振动模态参数的辨识值；

步骤4.2、判断振动模态参数的理论值和真实值之间的差别，若振动模态参数的理论值和辨识值之间相差大于10％，则根据辨识结果重新计算ZVD成型器的参数，并通过在轨编程方式修正存储在星载计算机中的ZVD成型器参数；如果理论值和辨识值相差小于等于10％，则不进行星上修正。

2.如权利要求1所述的主动抑制太阳帆板挠性振动的方法，其特征在于，所述的步进驱动电机采用正余弦细分驱动的步进电机。