飞轮振源特征参数提取方法及其分布特性

罗睿智，张激扬，樊亚洪，王晓伟，吴金涛

(1.北京控制工程研究所，北京 100094；2.中国空间技术研究院，北京 100094)

0 引 言

飞轮因具有结构紧凑、不消耗工质和寿命长等众多优点而被广泛应用于卫星等中小型航天器上，实现对航天器的三轴姿态机动与稳定控制。飞轮的内部结构[1]如图1所示，其中，轴承组件是飞轮的核心，其内部存在众多活动部件。在飞轮转子高速旋转的过程中，预紧轴承[2]的滚动面波纹[3]会引起预紧力波动[4]。它与飞轮转子不平衡产生的离心力和力偶共同构成飞轮微振动的两个主要振源因素，可能激励飞轮的模态和加剧转子涡动，从而引发结构共振。飞轮转子高速旋转时产生的宽频微振动通过飞轮基座传递至航天器舱板，成为影响航天器姿态稳定和指向精度的最大干扰源[5-6]。随着用户对航天器功能性能需求的提高，飞轮微振动对航天器平台稳定性的影响日益突显出来，制约着载荷性能的进一步提升，因此针对飞轮微振动的动力学特性研究[4]和抑制[7-8]成为惯性技术发展的重要方向。

文献[9]分析了转子振动加速度信号的频谱，指出飞轮径向振动频谱中主要包括工频及高次倍频成分。文献[10]等研究了飞轮的动力学模型，将飞轮的振源简化为离散的谐波。文献[11-12]拟合了飞轮的多阶谐波。可见对飞轮微振动研究主要集中于飞轮基座上输出的微振动频率成分的来源定位和幅值拟合。它们是振源经过飞轮结构传递的结果，掩盖了振源信息。可是飞轮振源幅值信息是飞轮产品微振动抑制和超静性能提升的基础，有必要对其进行深入研究。比较简便有效的方法是通过统计已有的飞轮微振动测试数据，提取飞轮振源幅值特征参数。

通过将飞轮固定在刚度较大的六分量微振动测量台上测试飞轮振动力和力矩(如图2所示)，忽略测试系统挠性对飞轮中、低频微振动幅值的影响。测得某50 Nms飞轮的微振动力/力矩瀑布图如图3所示。其中，横轴为频率，纵轴为转速，竖轴为振幅。

由图3可知，飞轮微振动输出包含了一系列放射状倍频振动，当它们与转子涡动曲线和模态曲线耦合时，可能产生较大结构共振，使得振源幅值在不同的转速和频率点处存在不同的放大。为了从这些振动数据中提取飞轮振源特征参数，并研究其分布特性，本文首先利用鲁棒回归从25个50Nms飞轮的微振动瀑布图中辨识出各飞轮各阶次(即：微振动频率对旋转频率之比)的振幅对转速二次幂系数，将其作为该飞轮振源幅值特征参数；其次，对多个飞轮的同一阶次特征参数进行Weibull分布拟合，并计算一些近似均匀的概率点的分位数，它们既宏观地展示了飞轮振源幅值的分布特性，也可将其作为飞轮振源幅值的分级标准。期望能为飞轮不平衡量和轴承预紧力波动量的评估和筛选提供参考。

1 基于鲁棒回归的特征参数提取

通过对飞轮微振动频率成分进行理论分析和测试统计，结果显示飞轮的微振动阶次[4]固定，其微振动特性的差异主要是振幅差异，则针对飞轮主要阶次的微振动幅值统计具有可行性。而输出的振幅是振源幅值经过飞轮内部结构传递作用的结果。从传递角度分析，未激发结构共振的低频振动在结构中的振动传递率近似为1；激发结构共振的中频段将被放大；而高频段存在一定的衰减，这表明：

1)对于未激发结构共振的频率成分的振幅主要取决于飞轮不平衡量和轴承的滚动面波纹等因素的激励作用，它反映了飞轮在该阶次的激励幅值。鉴于飞轮转子在高速旋转过程中产生的各阶次微振动幅值是转速的二次幂函数[10-13]，不妨假设第j阶振幅Fj(Ω)=βjΩ2。可将该系数βj作为该飞轮在该阶次振源幅值的一个特征参数，它较好地表征了飞轮内部的激励幅值随转速的变化规律。因此只需辨识出飞轮的主要阶次的特征参数，即可比较准确地反映出飞轮的微振动激励特征。

2)低频振动更真实反映了飞轮的微振动激励特征，因此在进行数据统计时应赋予其更大的权重，削弱其对共振处微振动数据的依赖。鉴于普通最小二乘法(OLS)易受离群值(共振点的振动数据)的杠杆作用使得对数据的拟合欠佳，拟合得到的模型的解释性和预测性相对较差，而M估计[14]、中位数回归和鲁棒回归[15-16]等方法都具有较好的鲁棒性，对离群值不敏感。因此本文将采用鲁棒回归提取飞轮在主要频率点处的特征参数。

对于线性归回，不妨假设待估参量βk是n维列向量，Xk是mk×n观测矩阵，Yk是mk维观测列向量，Ek是mk维残差列向量，且Yk=Xkβk+Ek，则第k次迭代的残差加权(权矩阵为Λk)平方和：

(1)

由此可解得估计系数为：

(2)

若式(2)中Λk=I，则变为普通最小二乘回归：

(3)

则观测向量Yk的OLS估计值为：

(4)

(5)

依据残差可计算出残差中位数Mk为：

(6)

式中：median为中位数函数。

误差项的标准离差估计sk：

(7)

则其标准化残差rk为：

(8)

其中，εk为前一次回归的残差；tc为调节常数。若将加权最小二乘的权重取为残差的函数(对残差较大的数据赋予较小的权重，反之亦然)，其变为鲁棒回归。常用的权函数有Bisquare核函数、Cauchy核函数和Welsch核函数等，其中，Bisquare核函数：

(9)

(10)

计算流程如图4所示。由此可回归得到飞轮各阶次振源幅值的二次幂系数βj(j∈{1,2,…,∞}),将其作为飞轮的振源特征参数。当飞轮结构不变时,系数βj也将不变，它是飞轮的不变特征量。

2 基于Weibull分布的特征参数分级

50 Nms飞轮产品的生产线和工艺过程固化，产品的研制质量严格受控，其微振动特性参数相对稳定，所积累的多个产品的微振动数据具有一定的统计规律性，故可通过统计的方法研究飞轮的振源特征参数分布特性。但微振动特性受到零件参数、装配间隙、轴承预紧力等众多因素的影响，因此又很难从理论上证明飞轮的微振动特征参数一定服从某种概率分布，故需进行试探性拟合。鉴于飞轮特征参数都大于零，根据直方图可见，其分布形式与广义Gamma分布或者Weibull分布[17]近似。该分布参数的取值可通过拟合得到，可将其作为该飞轮的分布特征参数。广义Gamma分布的概率密度函数：

(11)

式中:x>0，ρ>0为尺度参数；α>0为功率参数；λ>0为形状参数。Г(λ)为Gamma函数。当λ=1时，即可得到Weibull分布的概率密度和分布函数：

(12)

(13)

利用该Weibull分布来拟合各阶次微振动幅值特征参数，求出分布参数。则其逆函数为：

(14)

为了研究飞轮振源特征参数的分布特性，不妨以近似均匀的概率点p∈{0.2,0.4,0.6,0.8,0.95}作为其微振动分级概率，飞轮等级与分位数的对应关系如表1所示,再用式(14)求出表1的分位数。

表1 微振动特性分级表Table 1 Rating table of micro vibration characteristics

3 数据统计结果

飞轮主要频率阶次统计结果列于表2第2列。针对主要阶次在微振动瀑布图(缓慢升速时的振动频谱曲线随转速变化形成的三维图)中搜索振动幅值随转速变化的曲线，再利用图4求得该飞轮该阶次的二次特征参数βj。如某飞轮的工频和3倍频在升速过程中的振幅曲线及其拟合曲线如图5所示。

图5中曲线上振动较大的尖峰主要是因为该阶次振源激励被飞轮的结构共振放大所致。对未激发结构共振的工频振动的幅-频曲线具有较好的拟合精度，多次拟合曲线基本重合；对于存在结构共振的幅-频曲线，需采用图4的计算流程逐步消除结构共振导致的离群值，拟合曲线终将收敛于该阶次的激励幅值。该拟合曲线在一定程度上直观反映了该飞轮的微振动激励幅值随转速变化的关系。

利用图4的计算流程逐一提取飞轮各主要阶次的微振动激励幅值特征参数，进而可得到飞轮在全转速范围内的激励幅值-频率-转速三维图(某飞轮的激励三维图如图6所示)。与振动输出(图3所示)对比可见，回归得到的激励三维图消除了结构共振影响，直观反映了飞轮微振动激励幅值对频率和转速的关系。结果表明该飞轮除了工频(不平衡量导致的离心力/力偶)幅值较大以外，其余由轴承激发的预紧力波动量的幅值较小。可见，飞轮中、高频振动主要取决于这些微小振动被结构共振放大的倍数。

从25个50 Nms飞轮的微振动瀑布图数据中逐一提取各飞轮各主要阶次的振源激励特征参数，如表2所示。对表2中各阶次的所有特征参数按照Weibull分布进行拟合，从而得到各主要阶次的分布参数，再根据式(14)计算表1要求的分位数。可得到各主要阶次微振动激励幅值的分级图如图7所示，图中将飞轮的振动激励幅值分为六个等级，其中Ⅵ级为超出Ⅴ级以上的广泛区域。

从表2及图6和图7可以看出，飞轮的不平衡特征参数较大，且分布范围较宽，主要位于10-8量级；而飞轮轴承预紧力波动量的主要阶次的特征参数主要位于10-9量级。除{0.6,1,1.78,2.23,3.43,5.59}倍频的振动较大以外，其它倍频振动相对较小，但由于飞轮的结构阻尼也较小，这些微小的倍频振动也可能会激发结构共振，从而输出较大的振动。因此抑制飞轮的振动可主要从如下三个方面着手：1)减小飞轮的剩余不平衡量；2)增加飞轮的结构阻尼；3)研究飞轮的共振机理，优化模态曲线和涡动曲线分布，使之避开振幅较大的激励。

表2 25个飞轮在主要阶次的特征参数统计表Table 2 Characteristic parameters statistics at some main orders of 25 flywheels

4 结 论

为了研究飞轮振源幅值的分布特性，本文利用鲁棒回归的方法提取了与振源幅值直接相关的转速二次幂系数，将其作为飞轮的振源特征参数；并运用统计的方法对同阶次的特征参数进行拟合，计算出特定概率点的分位数。将该方法成功应用于对50 Nms飞轮的微振动数据统计分析中，从微振动瀑布图中提取了飞轮振源幅值特征参数，所计算出的分位数直观反映了该型飞轮的特征参数的分布情况。统计结果亦可作为该型飞轮振源特征参数的分级标准，为超静飞轮的筛选提供初步的判定准则，为飞轮改进提供参考依据，为飞轮在某种试验前后的结构一致性比对提供一种简单有效的方法。

[1] Sathyan K, Gopinath K, Hsu H Y, et al. Development of a lubrication system for momentum wheels used in spacecrafts [J]. Tribology Letters, 2008, 32(2):99-107.

[2] 王虹,卿涛,周刚,等.预载对飞轮用轴承组件振动特性的影响研究[J].空间控制技术与应用, 2015, 41(1):41-45. [Wang Hong, Qing Tao, Zhou Gang, et al. Effects of preload on the micro-vibration properties of bearing assemblies in flywheels [J]. Aerospace Control and Application, 2015, 41(1): 41-45.]

[3] Bai C Q, Xu Q Y. Dynamic model of ball bearings with internal clearance and waviness [J]. Journal of Sound and Vibration, 2006, 294(1-2):23-48.

[4] 罗睿智,虎刚,王全武.悬臂式SGCMG的高速转子的径向振动特性研究[J].振动与冲击, 2014(5):200-204. [Luo Rui-zhi, Hu Gang, Wang Quan-wu. Radial vibration characteristics of a high-speed rotor of a cantilevered SGCMG [J]. Journal of Vibration and Shock, 2014(5): 200-204.]

[5] 刘天雄,范本尧,杨慧.卫星飞轮扰振控制技术研究[J].航天器工程,2009,18(1):53-60. [Liu Tian-xiong, Fan Ben-yao, Yang Hui. Research on vibration disturbance control of momentum wheel [J]. Spacecraft Engineering, 2009, 18(1): 53-60.]

[6] 庞世伟,杨雷,曲广吉.高精度航天器微振动建模与评估技术最近进展[J].强度与环境,2007, 34(6):1-9. [Pang Shi-wei, Yang lei, Qu Guang-ji. New development of micro-vibration integrated modeling and assessment technology for high performance spacecraft[J]. Structure & Enviroment Enginee-ring, 2007, 34(6):1-9.]

[7] 刘天雄,林益明,王明宇,等.航天器振动控制技术进展[J]. 宇航学报, 2008, 29(1): 1-12. [Liu Tian-xiong, Lin Yi-ming, Wang Ming-yu, et al. Review of the spacecraft vibration control technology [J]. Journal of Astronautics, 2008, 29(1): 1-12.]

[8] 刘海平,王耀兵,史文华. 非线性能量阱对飞轮扰振特性的抑制[J]. 宇航学报, 2017, 38(5): 490-496. [Liu Hai-ping, Wang Yao-bing, Shi Wen-hua. Vibration suppression for a flywheel based on nonlinear energy sink [J]. Journal of Astronautics, 2017, 38(5): 490-496.]

[9] 王全武,虎刚.飞轮振动频谱特征的初步理论分析和验证[J].空间控制技术与应用, 2008, 34(4):42-46. [Wang Quan-wu, Hu Gang. Analysis and validation of flywheel vibration spectrum characteristics [J]. Aerospace Control and Application, 2008, 34(4):42-46.]

[10] 徐赵东,翁沉卉,朱俊涛.航天器中反作用轮干扰力仿真研究[J].振动、测试与诊断, 2013, 33(5):881-885. [Xu Zhao-dong, Weng Chen-hui, Zhu Jun-tao. Disturbance force simulation of reaction wheel in spacecraft [J]. Journal of Vibration, Measurement & Diagnosis, 2013, 33(5): 881-885.]

[11] Masterson R A, Miller D W, Grogan R L. Development and validation of reaction wheel disturbance models: empirical model [J]. Journal of Sound and Vibration, 2002, 249(3): 575-598.

[12] Masterson R A, Miller D W, Grogan R L. Development of empirical and analytical reaction wheel disturbance models [C]. The 40th Structures, Structural Dynamics, and Materials Conference and Exhibit, St. Louis, USA, April 12-15, 1999.

[13] Bialke B. High fidelity mathematical modeling of reaction wheel performance [C]. 1998 Annual AAS Rocky Mountain Guidance and Control Conference, Breckenridge, USA, Feb 4-8, 1998.

[14] Pennacchi P， Chatterton S， Ricci R. Rotor balancing using high breakdown-point and bounded-influence estimators [J]. Mechanical Systems and Signal Processing, 2010, 24(3): 860-872.

[15] 赵兰敏,郄志红,吴鑫淼,等.鲁棒回归在水工建筑物安全监控模型建立中的应用[J].中国农村水利水电,2007(3), 97-99,102. [Zhao Lan-min, Qie Zhi-hong, Wu Xin-miao, et al. Application of robust regression in the hydraulic structure safety monitoring modeling [J]. China Rural Water and Hydropower, 2007(3), 97-99,102.]

[16] 艾永乐,王玉栋,都静静,等. 基于LTS初值的稳健回归的谐波发射水平评估方法[J].电力系统保护与控制, 2015, 43(21):99-105. [Ai Yong-le, Wang Yu-dong, Du Jing-jing, et al. A method for assessing harmonic emission level based on robust regression of LTS initial value [J]. Power System Protection and Control, 2015, 43 (21):99-105.]

[17] 洪东跑,赵宇,马小兵. 广义Weibull分布参数的收缩估计[J].宇航学报, 2009, 30(6): 2442-2446. [Hong Dong-pao, Zhao yu, Ma Xiao-bing. Shrinkage estimation of generalized weibull parameter in censored samples [J]. Journal of Astronautics, 2009, 30(6): 2442-2446.]