-
摘要: 针对经验模态分解(Empirical mode decomposition, EMD)系列方法存在的模态分裂(Mode splitting, MS)问题, 提出中值互补集合经验模态分解(Median complementary ensemble EMD, MCEEMD)算法. 通过概率模型量化互补集合经验模态分解(Complementary ensemble EMD, CEEMD)的MS问题, 证明了使用中值算子替代算术平均算子对抑制MS的有效性. 为了兼具抑制MS和残留噪声的性能, MCEEMD算法首次在集合过程中结合了中值和平均算子. 具体地, 所提方法首先添加N对互补的白噪声至原信号中, 并经过EMD分解得到2N组固有模态函数(Intrinsic mode functions, IMFs), 然后分别对其中互补相关的IMFs两两取平均得到N组IMFs, 最后使用中值算子处理上述N组IMFs得到输出结果. 对仿真信号与两个真实案例的分析结果表明, 本文提出的MCEEMD方法不仅有效抑制了CEEMD的MS问题, 而且避免了单一使用中值算子的两个缺点: 分解完备性差和IMFs中存在的毛刺现象.
-
关键词:
- 模态分裂 /
- 中值算子 /
- 互补白噪声 /
- 互补集合经验模式分解
Abstract: In order to restrain the mode splitting (MS) problem of the empirical mode decomposition (EMD) series methods, a median complementary ensemble EMD (MCEEMD) algorithm is proposed in this paper. We first present novel probabilistic tools to quantify the MS phenomenon of complementary ensemble EMD (CEEMD), which aims at demonstrating the effectiveness of using the median operator to replace the mean operator during the ensemble process. To combine the advantages of suppressing MS and residual noise, the MCEEMD algorithm integrates both median and mean operators within the ensemble process for the first time. Specifically, the MCEEMD algorithm is enlightened and featured by following procedures: 1) Add N pairs of complementary white noise to the original signal to obtain 2N groups of intrinsic mode functions (IMFs) by EMD decomposition; 2) By averaging each pair of the complementary IMFs, the 2N groups of IMFs are computed into N IMF groups; 3) Assemble same-index components across the N groups of IMFs using the median operator to obtain the final IMFs within MCEEMD. Through typical simulations as well as two real-world cases, we show that the present work not only effectively alleviates the MS problem, but also avoids two shortcomings of using a single median operator, i.e., the poor decomposition completeness and the presence of burr in IMFs. -
随着航天技术的快速发展, 为满足不同需求, 人们开发了各种各样的航天器.其中, 姿轨控制系统是保障各类航天器稳定可靠运行的关键组成部分之一[1], 在较少的硬件配置和燃料消耗情况下, 能同时提高航天器的可靠性和降低航天器的综合成本, 这对实现航天器姿态的高效控制具有重要的意义.
航天器的姿态控制包含姿态速率阻尼和姿态机动, 主要通过喷气和飞轮控制两种方式实现.根据飞轮角动量交换控制原理, 一般需通过喷气控制或磁控制来阻尼航天器的转动, 而且喷气控制的硬件配置需超过4个.飞轮则适用于航天器姿态的精确控制, 通常会采用超过3个飞轮的硬件配置来应对飞轮长期运行中一些不确定因素引起的飞轮故障[2].当执行机构故障数目过多时, 航天器姿态控制会变为欠驱动系统[3-4], 此时航天器的姿态控制自由度少于姿态运动自由度, 导致无法使用常规的控制手段解决[5].即应用双飞轮的欠驱动姿态控制方案无法等同于全驱动控制方案, 当系统存在初始角动量时, 采用低于3个飞轮配置的航天器系统不可控[6].由此可见, 双飞轮控制方案存在非可控轴旋转时的角动量阻尼、动力学耦合效应对双飞轮控制的干扰和双飞轮无法直接实现非可控轴的姿态机动三个问题.另外, 虽然足够多的喷气配置能够实现航天器任意姿态下的姿态阻尼和姿态机动[7], 但这会增加航天器的燃料消耗, 并不利于降低航天器成本.
因此, 为了实现航天器姿态控制设计中的"低成本、高可靠"目标, 人们已经针对欠驱动航天器的燃料节省和姿态控制展开大量研究.段晨阳[7]针对喷气控制的配置与燃料节省问题, 提出了一种基于凸锥分析的过驱动航天器离线/在线混合推力动态分配法.郭朝礼等[8]针对具有推力器执行机构的航天器, 在一轴推力器失效的情况下, 设计了滑模变结构速率阻尼控制器.金磊等[9]提出了系统总角动量不为零时, 仅带两个飞轮的姿态稳定控制方法, 实现了航天器视线轴指向控制.戈新生等[10]根据航天器欠驱动轴的角速度分量是否为零分别设计相应的姿态稳定控制律, 提出了基于滑模函数的反馈控制律.而文献[11-13]主要分析了零角动量欠驱动卫星的控制问题, 分别针对单、双飞轮控制的卫星设计了分段光滑控制律.文献[5, 14]则以非零角动量欠驱动卫星为研究对象, 研究了不同稳定形式下的控制方法.宋道喆等[15]研究了轮控式零角动量欠驱动航天器姿态最优稳定控制问题.
本文为实现欠驱动角动量航天器的姿态控制, 提出共面双飞轮+单喷气的配置方案, 从欠驱动飞轮控制原理、喷气控制原理出发, 提出燃料最优角动量阻尼方法和全姿态机动方法.其中, 燃料最优角动量阻尼方法将航天器角动量消除, 全姿态机动方法在此基础上通过双飞轮实现三轴姿态机动.对具备一定角动量的航天器, 分析双飞轮姿态控制的稳定性, 介绍双飞轮实现喷气姿态阻尼燃料的最优方法和原理, 设计基于双轴旋转的双飞轮全姿态机动控制策略和控制律, 最后通过数学仿真验证控制方法的有效性.
1. 双飞轮-单喷气最优姿态阻尼原理
1.1 喷气最优组合控制
假设航天器装有$n$个喷气阀, 则喷气控制问题可转为对$n$个喷气阀的逻辑控制问题.其思路是将$n$个喷气阀的控制输出角动量与期望输出角动量等效, 最少燃料消耗的喷气阀组合为最优组合.
以$\Delta t_i$表示各个喷气阀的开机脉宽, 以$\bm{m}_i$表示各个喷气阀单位时间产生的角动量反向, 则航天器角动量阻尼控制问题可描述为
$$ \begin{equation} \sum\limits_{i=1}^{n}{\Delta {{t}_{i}}\cdot {{\bm{m}}_{i}}}=\bm{p} \end{equation} $$ (1) 式中, $\bm{p}$代表航天器本身的角动量.喷气阀的开机脉宽($\Delta t_i$)均非负, 至少分配4个喷气阀来实现全空间的角动量等效. $n\ (n\geq 4)$个喷气阀的开机脉宽调节可使得:
$$ \begin{equation} {f_{\min }} = {\rm{min}}\left( f \right), \quad f = \sum\limits_{i = 1}^n {\Delta {t_i}} \end{equation} $$ (2) 式中, $f$表示$n$个喷气阀的开机时间总和, $f_{\min}$表示诸多满足式(2)中的最小值, 该值越小, 表示喷气控制时所消耗的燃料越少.在至少4台推力器合理布局的前提下, 必定存在3个喷气阀(记为第$j$、$k$、$l$个)可以实现角动量等效:
$$ \begin{equation} \Delta {t_j}{{\bm{m}}_j} + \Delta {t_k}{{\bm{m}}_k} + \Delta {t_l}{{\bm{m}}_l} = {\bm{p}}, \quad \Delta t \ge 0 \end{equation} $$ (3) 当式(3)中向量$\bm{m}_j$、$\bm{m}_k$、$\bm{m}_l$线性无关时, 存在满足该等式的唯一解$\Delta t_j$、$\Delta t_k$、$\Delta t_l$, 必定存在一组解满足:
$$ \begin{equation} \left\{ \begin{array}{l} \Delta {t_j}{{\bm{m}}_j} + \Delta {t_k}{{\bm{m}}_k} + \Delta {t_l}{{\bm{m}}_l} = {\bm{p}}, \quad \Delta t \ge 0\\ {f_{\min }} = \Delta {t_j} + \Delta {t_k} + \Delta {t_l} \end{array} \right. \end{equation} $$ (4) 该式表明一定有对应第$j$、$k$、$l$个向量的解, 使得($\Delta t_j + \Delta t_k + \Delta t_l$)达到最小, 此最优问题可表示为:
$$ \begin{equation} \left\{ \begin{array}{l} \sum\limits_{i = 1}^n {\Delta {{\tilde t}_i}{{\bm{m}}_i}} = {\bm{p}}\\ \sum\limits_{i = 1}^n {\Delta {{\tilde t}_i}} \ge \left( {\Delta {t_j} + \Delta {t_k} + \Delta {t_l}} \right) = {f_{\min }} \end{array} \right. \end{equation} $$ (5) 式中, $\Delta {{\tilde{t}}_{i}}$表示第$i$个喷嘴的开机时长.由于$\bm{m}_j$、$\bm{m}_k$、$\bm{m}_l$在三维空间中线性无关, 对于任意的, 必然存在系数$\alpha_i$、$\beta_i$、$\gamma_i$满足:
$$ \begin{equation} {{\bm{m}}_{i}}={{\alpha }_{i}}{{\bm{m}}_{j}}+{{\beta }_{i}}{{\bm{m}}_{k}}+{{\gamma }_{i}}{{\bm{m}}_{l}} \end{equation} $$ (6) 结合式(5)和式(6)可得:
$$ \begin{equation} \sum\limits_{i=1}^{n}{\Delta {{{\tilde{t}}}_{i}}\left( {{\alpha }_{i}}{{\bm{m}}_{j}}+{{\beta }_{i}}{{\bm{m}}_{k}}+{{\gamma }_{i}}{{\bm{m}}_{l}} \right)}=\bm{p} \end{equation} $$ (7) 对比式(3)和式(7)可得:
$$ \begin{align} \sum\limits_{i=1}^{n}{\Delta {{{\tilde{t}}}_{i}}}=\, &\left( \Delta {{t}_{j}}+\Delta {{t}_{k}}+\Delta {{t}_{l}} \right)+ \nonumber\\ & \sum\limits_{i=1, i\ne j, k, l}^{n}{\Delta {{{\tilde{t}}}_{i}}\left( 1-{{\alpha }_{i}}-{{\beta }_{i}}-{{\gamma }_{i}} \right)} \end{align} $$ (8) 从式(8)可知, 当对于所有的$(1-\alpha_i-\beta_i-\gamma_i)$均大于或等于零时, 对于任意满足条件的时间总量$\sum\limits_{i=1}^{n}{\Delta {{{\tilde{t}}}_{i}}}$必定大于等于$(\Delta t_j + \Delta t_k + \Delta t_l)$, 即可证明该解是严格最优.下面可用反证法证明结论成立, 对于任意的$\bm{m}_i$, 必然存在以下等式关系:
$$ \begin{equation} \Delta {{t}_{j}}{{\bm{m}}_{j}}+\Delta {{t}_{k}}{{\bm{m}}_{k}}+\Delta {{t}_{l}}{{\bm{m}}_{l}}-\delta {{\bm{m}}_{i}}+\delta {{\bm{m}}_{i}}=\bm{p} \end{equation} $$ (9) 式中, $\delta$为正数且为小量.上式展开可写为:
$$ \begin{align} & \left( \vartriangle {{t}_{j}}-\delta {{\alpha }_{i}} \right) {{\bm{m}}_{j}}+\left( \vartriangle {{t}_{k}}-\delta {{\beta }_{i}} \right){{\bm{m}}_{k}} +\nonumber\\ &\qquad \left( \vartriangle {{t}_{l}}-\delta {{\gamma }_{i}} \right){{\bm{m}}_{l}}+\delta {{\bm{m}}_{i}}=\bm{p} \end{align} $$ (10) 若假设$(1-\alpha_i-\beta_i-\gamma_i)$小于零, 则$\alpha_i$、$\beta_i$、$\gamma_i)$中必有正数, 这就说明随着$\delta$增加, 上式中$\bm{m}_j$、$\bm{m}_k$、$\bm{m}_l$的系数会有一个先达到零.不妨令$\bm{m}_j$的系数先到零, 此时满足$\Delta t_j = \delta \alpha_i$, 即$\delta = \Delta t_j / \alpha_i$, 代入到式(10)可得:
$$ \begin{align} & \left( \Delta {{t}_{k}}-\frac{\Delta {{t}_{j}}}{{{\alpha }_{i}}}{{\beta }_{i}} \right){{\bm{m}}_{k}}+\nonumber \\ &\qquad \left( \Delta {{t}_{l}}-\frac{\Delta {{t}_{j}}}{{{\alpha }_{i}}}{{\gamma }_{i}} \right){{\bm{m}}_{l}}+\frac{\Delta {{t}_{j}}}{{{\alpha }_{i}}}{{\bm{m}}_{i}}=\bm{p} \end{align} $$ (11) 这意味着上式向量系数之和为:
$$ \begin{align} \Delta {{t}_{k}}&+\Delta {{t}_{l}}+\frac{\Delta {{t}_{j}}\left( 1-{{\beta }_{i}}-{{\gamma }_{i}} \right)} {{{\alpha }_{i}}} <\nonumber\\ & \Delta {{t}_{k}}+\Delta {{t}_{l}}+\Delta {{t}_{j}} \end{align} $$ (12) 这与$(\Delta t_j + \Delta t_k + \Delta t_l)$达到最小矛盾, 故对于所有的$(1-\alpha_i-\beta_i-\gamma_i)$均大于等于零.然而, 不受外力矩时航天器的角动量方向在惯性空间固定, 但航天器的姿态可变.如果航天器姿态能够按预期的姿态进行机动, 则能够产生诸多不同的喷气组合来实现最优阻尼.当航天器姿态机动至喷气阀中$\bm{m}_i$模值最大者的方向与航天器角动量方向重合时, 满足:
$$ \begin{equation} \Delta {{t}_{i}}{{\bm{m}}_{i}}=\bm{p} \end{equation} $$ (13) 式中, $\Delta t_i$表示喷气阀$i$的开机时长, 结合式(3)有:
$$ \begin{equation} \Delta {{t}_{i}}{{\bm{m}}_{i}}=\Delta {{t}_{j}}{{\bm{m}}_{j}}+\Delta {{t}_{k}}{{\bm{m}}_{k}}+\Delta {{t}_{l}}{{\bm{m}}_{l}} \end{equation} $$ (14) 由于$\bm{m}_i$模最大, 可得:
$$ \begin{align} \Delta {{t}_{i}}=\, &\Delta {{t}_{j}}\frac{\left( {{\bm{m}}_{j}}\cdot {{\bm{m}}_{i}} \right)}{\left( {{\bm{m}}_{i}}\cdot {{\bm{m}}_{i}} \right)}+\Delta {{t}_{k}}\frac{\left( {{\bm{m}}_{k}}\cdot {{\bm{m}}_{i}} \right)}{\left( {{\bm{m}}_{i}}\cdot {{\bm{m}}_{i}} \right)}+\nonumber \\ & \Delta {{t}_{l}}\frac{\left( {{\bm{m}}_{l}}\cdot {{\bm{m}}_{i}} \right)}{\left( {{\bm{m}}_{i}}\cdot {{\bm{m}}_{i}} \right)}\le \Delta {{t}_{j}}+\Delta {{t}_{k}}+\Delta {{t}_{l}} \end{align} $$ (15) 上式等号当且仅当$\bm{m}_j = \bm{m}_k = \bm{m}_i = \bm{m}_i$时成立, 这表明只需将航天器的$\bm{m}_i$机动至与航天器角动量方向重合便可实现燃料消耗的全局最优.
因此, 对于给定姿态, 一定存在3个喷气阀的组合方案实现燃料消耗最优.而对于非给定姿态, 当航天器角动量与喷气阀作用力矩反向同向时, 便可实现全局燃料最优, 此时双飞轮组合能够控制喷气作用力矩反向与航天器角动量同向.
1.2 双飞轮控制原理
双飞轮仅能直接产生两个方向的控制力矩, 在航天器转动时还会伴随耦合力矩, 该耦合力矩可用来阻尼航天器的角速度.
双飞轮进行航天器的角动量阻尼控制的方法为:
$$ \begin{equation} \left( \begin{aligned} & {{{\dot{h}}}_{1}} \\ & {{{\dot{h}}}_{2}} \\ \end{aligned} \right)=\left( \begin{aligned} & {{d}_{1}}\bm{\omega }\cdot {{\bm{e}}_{1}} \\ & {{d}_{2}}\bm{\omega }\cdot {{\bm{e}}_{2}} \\ \end{aligned} \right), \quad {{d}_{1}}>0, {{d}_{2}}>0 \end{equation} $$ (16) 式中, ${{\dot{h}}_{1}}$和${{\dot{h}}_{2}}$分别表示飞轮1和2的角动量变化率, $\omega$表示航天器的角速度. $\bm{e}_1$和$\bm{e}_2$分别表示飞轮1和2的旋转方向. $d_1$、$d_2$为控制系数.上式表明飞轮轴向的航天器角速度会被阻尼, 而在双飞轮无法控制的方向上, 航天器的角速度不会直接受飞轮阻尼而消除.航天器的动力学方程可以写为:
$$ \begin{equation} {J}{\dot{\bm\omega }}+\bm{\omega }\times \left( {J}\bm{\omega }+{{\bm{h}}_{1}}+{{\bm{h}}_{2}} \right)+{{{\dot{\bm h}}}_{1}}+{{{\dot{\bm h}}}_{2}}=0 \end{equation} $$ (17) 将双飞轮分别安装在航天器的$x$、$y$轴上, 上式写成分量形式:
$$ \begin{equation} \begin{cases} {{J}_{11}}{{{\dot{\omega }}}_{x}}+\left( {{J}_{33}}-{{J}_{22}} \right){{\omega }_{y}}{{\omega }_{z}}+{{d}_{1}}{{\omega }_{x}}- {{h}_{2}}{{\omega }_{z}}=0 \\ {{J}_{22}}{{{\dot{\omega }}}_{y}}+\left( {{J}_{11}}-{{J}_{33}} \right){{\omega }_{x}}{{\omega }_{z}}+{{d}_{2}}{{\omega }_{y}}+ {{h}_{1}}{{\omega }_{z}}=0 \\ {{J}_{33}}{{{\dot{\omega }}}_{z}}+\left( {{J}_{22}}-{{J}_{11}} \right){{\omega }_{x}}{{\omega }_{y}}+{{h}_{2}}{{\omega }_{x}}- {{h}_{1}}{{\omega }_{y}}=0 \end{cases} \end{equation} $$ (18) 对上式左右点乘角速度矢量$\bm{\omega}$有:
$$ \begin{align} {{J}_{11}}{{\omega }_{x}}{{{\dot{\omega }}}_{x}}&+{{J}_{22}}{{\omega }_{y}}{{{\dot{\omega }}}_{y}}+{{J}_{33}}{{\omega }_{z}}{{{\dot{\omega }}}_{z}}= \nonumber\\ & -{{d}_{1}}{\omega }_{x}^{2}-{{d}_{2}}{\omega }_{y}^{2} \end{align} $$ (19) 以$\Psi =\frac{1}{2}\left({{J}_{11}}{\omega }_{x}^{2}+{{J}_{22}}{\omega }_{y}^{2}+{{J}_{33}}{\omega }_{z}^{2} \right)$表示航天器自身角速度的能量函数, 上式可以写成:
$$ \begin{equation} \dot{\Psi }=-{{d}_{1}}{\omega }_{x}^{2}-{{d}_{2}}{\omega }_{y}^{2}\le 0 \end{equation} $$ (20) 此处定义航天器维持仅$z$轴方向的旋转且飞轮角动量均为零的状态为奇异状态, 奇异状态不受外部干扰时会维持.下面说明三个问题:
1) 奇异状态是一个不稳定状态
假设对奇异状态施加飞轮小转动的扰动$\Delta h$时, 航天器、飞轮组成的系统不再是奇异状态, 航天器的能量函数$\Psi $会衰减, 用$\omega_{z0}$表示奇异状态的航天器$z$轴角速度, 则:
$$ \begin{equation} \frac{1}{2}\left( {{J}_{11}}{\omega }_{x}^{2}+{{J}_{22}}{\omega }_{y}^{2}+{{J}_{33}}{\omega }_{z}^{2} \right)<\frac{1}{2}\left( {{J}_{33}}{\omega }_{z0}^{2} \right) \end{equation} $$ (21) 假设航天器还会回到奇异状态, 则由角动量守恒定律:
$$ \begin{equation} {{J}_{33}}{{\omega }_{z}}=\sqrt{{{\left( {{J}_{33}}{{\omega }_{z0}} \right)}^{2}}+\Delta {{h}^{2}}} \end{equation} $$ (22) 很明显式(21)与式(22)矛盾.所以, 奇异状态是一个不稳定状态.
2) 航天器停转状态是一个稳定状态
航天器停转时, 航天器的角动量都被飞轮吸收, 航天器的角速度能量函数降到最低0.当航天器不在停转状态和奇异状态时, 分为两种情况分析稳定性.
第一种情况为${\omega}_z = 0$时, 因航天器不在停转状态, 所以${\omega}_x$、${\omega}_y$不会同时为零, 即能量函数的导数小于零, 能量函数会衰减.
第二种情况为${\omega}_z \ne 0$时, 因航天器不在奇异状态, 所以${\omega}_x = {\omega}_y = {h}_1 = {h}_2=0$不会成立.只要${\omega}_x$、${\omega}_y$中至少有一个非零, 则能量函数会衰减; 在${\omega}_x$、${\omega}_y$均为零时, 则${h}_1$、${h}_2$中至少有一个非零, 从式(18)中看出, ${\omega}_x$、${\omega}_y$的导数至少有一个非零, 这就意味着${\omega}_x$、${\omega}_y$均为零的状态会立即消失, 能量函数最终还是会衰减.
因此航天器不在停转状态和奇异状态时, 其能量函数会衰减, 同时航天器的奇异状态是一个不稳定状态, 因此航天器停转状态是一个全局稳定状态.
3) 阻尼系数较大、转动惯量均匀时, $\omega_z$几乎单调衰减
在阻尼系数较大、转动惯量均匀时, 航天器的$\omega_x$、$\omega_y$会被快速阻尼至零, 而其外激励会使之发生变化.由式(18), 认为转动惯量均匀, 近似有:
$$ \begin{equation} \left\{ \begin{aligned} & {{\omega }_{x}}\approx \frac{1}{{{J}_{11}}{s}+{{d}_{1}}}\left( {{h}_{2}}{{\omega }_{z}} \right) \\ & {{\omega }_{y}}\approx \frac{1}{{{J}_{22}}{s}+{{d}_{2}}}\left( -{{h}_{1}}{{\omega }_{z}} \right) \\ \end{aligned} \right. \end{equation} $$ (23) 因此, 可以认为$\omega_x$和$\omega_y$分别是$h_2 \omega_z$和$-h_1 \omega_z$延迟响应, 即它们同号.更进一步根据${{J}_{33}}{{\dot{\omega }}_{z}}\approx {{h}_{1}}{{\omega }_{y}}-{{h}_{2}}{{\omega }_{x}}\approx -{{\lambda }_{1}}{{h}_{1}}^{2}{{\omega }_{z}}-{{\lambda }_{2}}{{h}_{2}}^{2}{{\omega }_{z}}$可知$\omega_z$是单调衰减的.
综上所述, 只要系统不属于仅绕$z$轴旋转且飞轮不转的奇异状态, 则最终$\omega_x$、$\omega_y$、$\omega_z$趋于零, 系统角动量转移到了双飞轮角动量上.此时, 可通过安装在该平面上的单喷气作用力矩反向机动至航天器角动量方向进行喷气阻尼.由于航天器仅在$xoy$平面上有执行能力, 应该通过绕$xoy$平面上某个轴进行转动, 如图 1所示.
以$\bm{H} = (h_x, h_y, h_z)$表示航天器角动量方向, $\bm{m} = (m_x, m_y, m_z)$表示喷气在本体下的作用力矩反方向, 不妨令$\bm{m}$为$x$轴, 则在$xoy$平面上的机动轴$\bm{q} = (q_x, q_y, q_z)$可以表示为:
$$ \begin{equation} \begin{cases} {{q}_{x}}=\frac{{{h}_{y}}}{{{h}_{x}}^{2}+{{h}_{y}}^{2}+{{h}_{z}}^{2}} \\ {{q}_{y}}=1-\frac{{{h}_{x}}}{{{h}_{x}}^{2}+{{h}_{y}}^{2}+{{h}_{z}}^{2}} \\ {{q}_{z}}=0 \end{cases} \end{equation} $$ (24) 控制律的设计为:
$$ \begin{equation} \left\{ \begin{aligned} & \kappa =\sqrt{{{\left( {{h}_{z}}{{q}_{y}} \right)}^{2}}+{{\left( {{h}_{z}}{{q}_{x}} \right)}^{2}}+{{\left( {{h}_{y}}{{q}_{x}}-{{h}_{x}}{{q}_{y}} \right)}^{2}}} \\ & {{T}_{x}}=-{{D}_{x}}{{\omega }_{x}}-{{P}_{x}}\frac{{{h}_{z}}{{q}_{x}}{{q}_{y}}}{\left| {{q}_{y}} \right|\kappa } \\ & {{T}_{y}}=-{{D}_{y}}{{\omega }_{y}}-{{P}_{y}}\frac{{{h}_{z}}{{q}_{y}}^{2}}{\left| {{q}_{y}} \right|\kappa } \\ \end{aligned} \right. \end{equation} $$ (25) 式中, $T_x$、$T_y$为控制力矩, $D_x$、$D_y$、$P_x$、$P_y$为控制系数.喷气力矩反向与航天器角动量方向接近时, 使用喷气进行阻尼, 控制律为:
$$ \begin{equation} {{\bm{T}}_{jet}}=-\xi \bm{H} \end{equation} $$ (26) 式中, $\bm{T}_{jet}$为喷气控制力矩, $\xi$为控制系数.
基于喷气最优组合和双飞轮控制原理, 通过理论推导给出共面双飞轮-单喷气的航天器姿态阻尼控制方案, 可在消耗燃料最少时, 实现航天器的姿态阻尼控制.当姿态阻尼完成后, 会得到的零角动量航天器, 进而使用双飞轮进行的全姿态机动控制.
2. 全姿态机动原理
第1节描述了双飞轮单轴姿态机动及喷气最优阻尼原理, 本节将在姿态阻尼完成之后得到零角动量航天器, 以双飞轮完成全姿态机动(三轴姿态机动)的控制设计.
姿态机动的本质可以描述为航天器上的两个固连矢量指向预期方位.可通过$xoy$平面上的转动轴实现固连矢量指向控制, 分为绕第一个轴转动和绕第二个轴转动.在第一个固连矢量机动至与预期方位重合时, 再绕第一个固连矢量旋转使第二个固连矢量机动到位, 如图 2所示:
基于双轴旋转的的控制方法, 可实现姿态机动, 用$\bm{e}_x=(1, 0, 0)$、$\bm{e}_y=(0, 1, 0)$表示航天器固连$x$轴、$y$轴, 用$\bm{e}_{xf}=(e_{xfx}, e_{xfy}, e_{xfz})$、$\bm{e}_{yf}=(e_{yf x}, e_{yf y}, e_{yf z})$表示对应$\bm{e}_x$、$\bm{e}_y$的期望方向:
$$ \begin{equation} \left\{ \begin{aligned} & {{\bm{C}}_{x}}={{\bm{e}}_{x}}\times {{\bm{e}}_{xf}} \\ & {{\bm{C}}_{y}}={{\bm{e}}_{y}}\times {{\bm{e}}_{yf}} \\ \end{aligned} \right. \end{equation} $$ (27) 如果$\bm{C}_x$更接近$xoy$平面, 则$\bm{q}=(q_x, q_y, q_z)$为:
$$ \begin{equation} \left\{ \begin{aligned} & {{q}_{x}}={{e}_{xfy}} \\ & {{q}_{y}}=1-{{e}_{xfx}} \\ & {{q}_{z}}=0 \\ \end{aligned} \right. \end{equation} $$ (28) $$ \begin{equation} \left\{ \begin{aligned} & \kappa =\sqrt{{{\left( {{e}_{xfz}}{{q}_{y}} \right)}^{2}}+{{\left( {{e}_{xfz}}{{q}_{x}} \right)}^{2}}+{{\left( {{e}_{xfy}}{{q}_{x}}-{{e}_{xfx}}{{q}_{y}} \right)}^{2}}} \\ & {{T}_{x}}=-{{D}_{x}}{{\omega }_{x}}-{{P}_{x}}\frac{{{e}_{xfz}}{{q}_{x}}{{q}_{y}}}{\left| {{q}_{y}} \right|\kappa } \\ & {{T}_{y}}=-{{D}_{y}}{{\omega }_{y}}-{{P}_{y}}\frac{{{e}_{xfz}}{{q}_{y}}^{2}}{\left| {{q}_{y}} \right|\kappa } \\ \end{aligned} \right. \end{equation} $$ (29) 如果$\bm{C}_y$更接近$xoy$平面, 则$\bm{q}=(q_x, q_y, q_z)$为:
$$ \begin{equation} \left\{ \begin{aligned} & {{q}_{x}}=1-{{e}_{yfy}} \\ & {{q}_{y}}={{e}_{yfx}} \\ & {{q}_{z}}=0 \\ \end{aligned} \right. \end{equation} $$ (30) $$ \begin{equation} \left\{ \begin{aligned} & \kappa =\sqrt{{{\left( {{e}_{yfz}}{{q}_{y}} \right)}^{2}}+{{\left( {{e}_{yfz}}{{q}_{x}} \right)}^{2}}+{{\left( {{e}_{yfy}}{{q}_{x}}-{{e}_{yfx}}{{q}_{y}} \right)}^{2}}} \\ & {{T}_{x}}=-{{D}_{x}}{{\omega }_{x}}+{{P}_{x}}\frac{{{e}_{yfz}}{{q}_{x}}^{2}}{\left| {{q}_{x}} \right|\kappa } \\ & {{T}_{y}}=-{{D}_{y}}{{\omega }_{y}}+{{P}_{y}}\frac{{{e}_{yfz}}{{q}_{x}}{{q}_{y}}}{\left| {{q}_{x}} \right|\kappa } \\ \end{aligned} \right. \end{equation} $$ (31) 基于双飞轮全姿态机动原理, 在两个固连矢量中, 通过分析它们机动至目标指向的转动角度大小, 可先使机动转动角度较小的固连矢量与预期方位接近重合, 再实现另一个固连矢量机动至预期方位.
3. 数值仿真验证及分析
为了验证上述提出的共平面双飞轮-单喷气最优控制方法, 仿真过程包括星体角动量阻尼、全姿态机动过程.为体现阻尼方法能够实现燃料最优, 对比力矩方向分别在$+x$、$+y$、$+z$、$-x$、$-y$、$-z$的6喷气组合, 以相同的初始条件进行仿真验证.本文以某航天器1°/s的初始速率、飞轮存在初始1 Nms的角动量进行仿真计算, 并将喷气安装在$+y$轴上, 使之能够产生$-x$方向的力矩, 如图 3所示.
1) 星体角动量阻尼仿真
针对双飞轮航天器阻尼, 理论分析发现可将航天器角速度阻尼至零附近.图 4给出了仿真中航天器角速度的变化示意图.发现通过施加双飞轮阻尼航天器转动的控制律, 可实现对$x$、$y$轴角速度的直接阻尼控制, 以及对$z$轴角速度的间接阻尼控制.从图中可以看出, 在控制律的作用下, $z$轴角速度呈现衰减的趋势, 最终航天器阻尼的精度达到了$4.0110 \times 10^{5}$ rad/s.航天器的角动量最终被飞轮吸收, 双飞轮最终将航天器$z$轴的角动量衰减至约$4 \times 10 ^5$ Nms.
图 4 双飞轮阻尼时航天器角速度的变化Fig. 4 The variation of spacecraft angular velocity by two flywheels' damping图 5给出了针对双飞轮姿态机动至喷气阻尼的仿真结果, 图中显示了双飞轮控制喷气力矩反方向至航天器角动量方向的过程及喷气在二者方向接近重合时进行阻尼控制的过程.控制过程中$x$轴角速度要大于其他两个轴, 当喷气阻尼动作执行后, 航天器的角速度得到衰减并趋于稳定.图 6为航天器总角动量的变化情况, 发现在喷气作用下, 航天器的角动量呈现阶梯式衰减, 最终衰减至$1$ $\times$ $10^3$ Nms达到稳定.
图 5 双飞轮姿态机动喷气阻尼的航天器角速度变化Fig. 5 The variation of spacecraft angular velocity by two flywheels' maneuver and jet's damping
图 6 双飞轮姿态机动喷气阻尼的航天器角动量变化Fig. 6 The variation of spacecraft angular momentum by two flywheels' maneuver and jet's damping图 7给出了双飞轮姿态机动喷气启动阻尼的工作时间和6喷气组合启动阻尼的工作时间变化的对比示意图, 当航天器稳定工作时, 航天器系统角动量消耗量为1.33 Nms, 喷气消耗的角动量为1.33 Nms, 二者相同, 表明喷气消耗的燃料无浪费.而使用6喷气组合直接进行姿态阻尼的喷气角动量消耗量为1.97 Nms, 增加了48.1 %, 可见双飞轮喷气阻尼的方案是最优的.
2) 全姿态机动仿真
当双飞轮姿态机动偏航90°时, 将阻尼控制及姿态机动的全过程进行仿真演示, 姿态角速度控制精度为0.1°/s, 仿真结果如图 8和9所示.稳定之后, 姿态接近目标值后, 机动轴在$x$轴与$y$轴间切换, 呈现出一定的震荡特性, 如图 8所示, 从图 9发现最终的姿态控制精度小于1°.
4. 结束语
本文针对欠驱动航天器的姿态控制, 结合喷气姿态阻尼的角动量等效原理, 在给定姿态、非给定姿态两种情况下推导出燃料消耗最小时喷气阀的最优组合方案.进一步为了实现喷气的全局最优, 结合欠驱动飞轮姿态控制策略, 分析欠驱动飞轮航天器的姿态控制原理及稳定性, 提出共面双飞轮-单喷气的配置方案, 该方案可通过双飞轮组合稳定航天器的角速度, 并达到燃料全局最省的预期姿态机动.在此基础上, 提出绕两轴旋转的姿态机动路径规划方法来实现航天器姿态机动与稳定控制.
最后, 对双飞轮-单喷气的姿态速率阻尼和姿态机动进行仿真, 验证了双飞轮-单喷气的姿态阻尼方法能够实现喷气燃料消耗最少; 验证了两次旋转轴机动的双飞轮控制方法能够实现姿态机动控制.仿真结果表明共面双飞轮-单喷气的欠驱动姿态阻尼方案能够在较少硬件配置下实现对航天器的姿态控制, 且消耗燃料最少.
另外, 本文研究中尚未考虑挠性航天器的欠驱动最优控制问题, 未来将在挠性振动干扰之下研究欠驱动双飞轮-单喷气航天器的最优控制特性.
-
图 1 EMD分解噪声辅助信号得到的前5个互补IMFs
Fig. 1 The first five complementary IMFs obtained from the noise-assisted signal through EMD
图 2 互补IMFs (由噪声辅助信号得到)的${p_i}(f)$和${r_i}(f)$曲线
Fig. 2 The curves ${p_i}(f)$ and ${r_i}(f)$ corresponding to the complementary IMFs (obtained from the noise-assisted signal)
图 3 互补IMFs (由噪声辅助信号得到)的${P_i}(f)$曲线
Fig. 3 The curves ${P_i}(f)$ corresponding to the complementary IMFs (obtained from the noise-assisted signal)
图 4 不同算子处理互补IMFs集合得到的$MSD(f)$曲线
Fig. 4 Curves of $MSD(f)$ obtained by processing the complementary IMFs with different operators
图 6 不同集合尺寸下CEEMD、MEEMD和MCEEMD的$MSR(f)$曲线
Fig. 6 $MSR(f)$curves for different ensemble sizes within CEEMD, MEEMD and MCEEMD
图 7 MCEEMD、MEEMD和CEEMD在不同集合尺寸下的$SDR(N)$曲线
Fig. 7 $SDR(N)$ curves for different ensemble sizes within MCEEMD, MEEMD and CEEMD
图 8 4种方法分解仿真信号所得的前5个IMF
Fig. 8 The first five IMFs obtained by decomposing the simulated signal by four methods
图 9 MEEMD、MCEEMD分解结果中的${d_{\rm{2}}}$模态
Fig. 9 The${d_{\rm{2}}}$mode in the decomposition results of MEEMD and MCEEMD
图 10 4种方法分解结果的PSD
Fig. 10 The PSD curves of the decomposition results from the four methods
图 11 MCEEMD分解血流信号所得的前8个分量
Fig. 11 The first eight components of the blood flow signal decomposed by MCEEMD
图 13 EEMD、CEEMD、MEEMD和MCEEMD提取的血流成分频率归一化功率谱
Fig. 13 Frequency normalized PSD of the blood flow component extracted by EEMD, CEEMD, MEEMD, and MCEEMD, respectively
表 1 4种方法的性能指标
Table 1 Performance indicators of the four methods
方法 PCC RMSE PSD面积比 (%) EEMD 0.9568 0.3087 0.28 CEEMD 0.9986 0.0031 0.21 MEEMD 0.7293 1.2938 0.24 MCEEMD 0.9980 0.1614 0.14 
下载: 导出CSV
表 2 4种方法的计算时间
Table 2 Calculation time of the four methods
方法 计算时间 (s) EEMD 14.32 CEEMD 28.95 MEEMD 14.58 MCEEMD 29.01
下载: 导出CSV
-
[1] Huang N E, Shen Z, Long S R, Wu M C, Shih H H, Zheng Q, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proceedings Mathematical Physical & Engineering ences, 1998, 454(1971): 903-995. doi: 10.1098/rspa.1998.0193 [2] Mandic D P, Rehman N U, Wu Z, Huang N E. Empirical mode decomposition-based time-frequency analysis of multivariate signals: the power of adaptive data analysis. Signal Processing Magazine IEEE, 2013, 30(6): 74-86. doi: 10.1109/MSP.2013.2267931 [3] 杨默涵, 陈万忠, 李明阳. 基于总体经验模态分解的多类特征的运动想象脑电识别方法研究. 自动化学报, 2017, 43(5): 743−752.Yang Mo-Han, Chen Wan-Zhong, Li Ming-Yang. Multiple feature extraction based on ensemble empirical mode decomposition for motor imagery EEG recognition tasks. Acta Automatica Sinica, 2017, 43(5): 743-752. [4] Lang X, Lu S, Xie L, Zakharov A, Zhong D, Sirkka-Liisa Jämsä-Jounela. Bihocerence based industrial control loop nonlinearity detection and diagnosis in short nonstationary time series. Journal of Process Control, 2018, 63: 15-28. doi: 10.1016/j.jprocont.2018.01.001 [5] Lang X, Zhang Y F, Xie L, Jin X, Horch A, Su H Y. Use of fast multivariate empirical mode decomposition for oscillation monitoring in noisy process plant. Industrial & Engineering Chemistry Research, 2020, 59(25): 11537-51. doi: 10.1021/acs.iecr.9b06351 [6] Devi A S, Maragatham G, Boopathi K, Rangaraj A G. Hourly day-ahead wind power forecasting with the EEMD-CSO-LSTM-EFG deep learning technique. Soft Computing, 2020, 24(16): 12391-12411. doi: 10.1007/s00500-020-04680-7 [7] 李霞, 卢官明, 闫静杰, 张正言. 多模态维度情感预测综述. 自动化学报, 2018, 44(12): 2142-2159.Li Xia, Lu Guan-Ming, Yan Jing-Jie, Zhang Zheng-Yan. A survey of dimensional emotion prediction by multimodal cues. Acta Automatica Sinica, 2018, 44(12): 2142-2159. [8] Wu Z, Huang N E. Ensemble empirical mode decomposition: a noise-assisted data analysis method. Advances in Adaptive Data Analysis, 2009, 1(1): 1-41. doi: 10.1142/S1793536909000047 [9] Flandrin P, Rilling G, Goncalves P. Empirical mode decomposition as a filter bank. IEEE signal processing letters, 2004, 11(2): 112-114. doi: 10.1109/LSP.2003.821662 [10] Yeh J R, Shieh J S, Huang N E. Complementary ensemble empirical mode decomposition: A novel noise enhanced data analysis method. Advances in Adaptive Data Analysis, 2010, 02(02): 135-156. doi: 10.1142/S1793536910000422 [11] Torres M E, Colominas M A, Schlotthauer G, Flandrin P. A complete ensemble empirical mode decomposition with adaptive noise. In: Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). Prague, Czech Republic: IEEE, 2011. 4144−4147 [12] Rehman N, Mandic D P. Multivariate empirical mode decomposition. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2010, 466(2117): 1291-1302. doi: 10.1098/rspa.2009.0502 [13] Lang X, Zheng Q, Zhang Z M, Lu S, Xie L. Fast multivariate empirical mode decomposition. IEEE Access, 2018, 6: 65521-65538. doi: 10.1109/ACCESS.2018.2877150 [14] Lang X, ur Rehman N, Zhang Y F, Xie L, Su H Y. Median Ensemble Empirical Mode Decomposition. Signal Processing, 2020, 176: 107686. doi: 10.1016/j.sigpro.2020.107686 [15] Cheng J S, Yu D J, Tang J S, Yang Y. Application of SVM and SVD technique based on EMD to the fault diagnosis of the rotating machinery. Shock and Vibration, 2009, 16(1): 89-98. doi: 10.1155/2009/519502 [16] 赵春晖, 余万科, 高福荣. 非平稳间歇过程数据解析与状态监控 —回顾与展望. 自动化学报, 2020, 46(10): 2072−2091. doi: 10.16383/j.aas.c190586Zhao Chun-Hui, Yu Wan-Ke, Gao Fu-Rong. Data analytics and condition monitoring methods for nonstationary batch processes — current status and future. Acta Automatica Sinica, 2020, 46(10): 2072−2091. doi: 10.16383/j.aas.c190586 [17] 周成江, 吴建德, 杨静宗. 基于CEEMD-SVD-LSSVM的矿浆管线核心设备故障诊断. 云南大学学报自然科学版, 2018, 40(05): 886-896.Zhou Cheng-Jiang, Wu Jian-De, Yang Jing-Zong. Fault diagnosis of core equipment of slurry pipeline based on CEEMD-SVD-LSSVM. Journal of Yunnan University Natural Sciences Edition, 2018, 40(05): 886-896. [18] 罗继辉, 黄国勇. 基于广义S变换和深度置信网络的单向阀故障诊断. 电子测量与仪器学报, 2019, 33 (09): 197-203.Luo Ji-Zhong, Huang Guo-Yong. Check valve fault diagnosis based on generalized S-transform and deep belief network. Journal of Electronic Measurement and Instrumentation, 2019, 33 (09): 197-203. [19] He B B, Zhang Y F, Zhang K X, Chen J H, Zhang J H, Liang H. Optimum speckle tracking based on ultrafast ultrasound for improving blood flow velocimetry. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2020, 68(3): 494-509. [20] 夏平, 施宇, 雷帮军, 龚国强, 胡蓉, 师冬霞. 复小波域混合概率图模型的超声医学图像分割. 自动化学报, 2021, 47(1): 185-196.Xia Ping, Shi Yu, Lei Bang-Jun, Gong Guo-Qiang, Hu Rong, Shi Dong-Xia. Ultrasound medical image segmentation based on hybrid probabilistic graphical model in complex-wavelet domain. Acta Automatica Sinica, 2021, 47(1): 185-196. [21] 林文晶, 张榆锋, 章克信, 李支尧, 李海燕, 高莲, 等. 总体经验模态细分法提取血流超声多普勒信号的研究. 电子学报, 2014, 42(007): 1424-1428.Lin Wen-Jing, Zhang Yu-Feng, Zhang Ke-Xin, Li Zhi-Yao, Li Hai-Yan, Gao Lian, et al. Extraction of doppler ultrasound blood signals using the delicate separation method based on the EEMD algorithm. Acta Electronica Sinica, 2014, 42(007): 1424-1428. [22] Zhang Y F, Gao Y L, Wang L, Chen J H, Shi X L. The removal of wall components in Doppler ultrasound signals by using the empirical mode decomposition algorithm. IEEE Transactions on Biomedical Engineering, 2007, 54(9): 1631-1642. doi: 10.1109/TBME.2007.891936 [23] Majd S M M T, Asl B M. Adaptive spectral doppler estimation based on the modified amplitude spectrum capon. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2021, 68(5): 1664-1675. doi: 10.1109/TUFFC.2020.3044774 [24] Lang X, Zheng Q, Xie L, Horch A, Su H Y. Direct Multivariate Intrinsic Time-Scale Decomposition for Oscillation Monitoring. IEEE Transactions on Control Systems Technology, 2020, 28(6): 2608-2615. doi: 10.1109/TCST.2019.2940374 [25] Xie X Y, Liu H, Shu M L, Zhu Q, Huang A P, Kong X Q, et al. A multi-stage denoising framework for ambulatory ECG signal based on domain knowledge and motion artifact detection. Future Generation Computer Systems, 2021, 116: 103-116. doi: 10.1016/j.future.2020.10.024 期刊类型引用(2)
1. 陈启明，文青松，郎恂，谢磊，苏宏业. 一元及多元信号分解发展历程与展望. 自动化学报. 2024(01): 1-20 . 
本站查看2. 郎煜民，郎恂，吴建德，刘燕，李鹏. 一种新的工业过程振荡数据去噪技术. 仪器仪表学报. 2024(12): 307-318 . 百度学术其他类型引用(2)
-
下载:
计量
- 文章访问数: 1842
- HTML全文浏览量: 556
- PDF下载量: 159
- 被引次数: 4
