飛輪振源特征參數(shù)提取方法及其分布特性

羅睿智，張激揚(yáng)，樊亞洪，王曉偉，吳金濤

(1.北京控制工程研究所，北京 100094；2.中國空間技術(shù)研究院，北京 100094)

0 引 言

飛輪因具有結(jié)構(gòu)緊湊、不消耗工質(zhì)和壽命長等眾多優(yōu)點而被廣泛應(yīng)用于衛(wèi)星等中小型航天器上，實現(xiàn)對航天器的三軸姿態(tài)機(jī)動與穩(wěn)定控制。飛輪的內(nèi)部結(jié)構(gòu)[1]如圖1所示，其中，軸承組件是飛輪的核心，其內(nèi)部存在眾多活動部件。在飛輪轉(zhuǎn)子高速旋轉(zhuǎn)的過程中，預(yù)緊軸承[2]的滾動面波紋[3]會引起預(yù)緊力波動[4]。它與飛輪轉(zhuǎn)子不平衡產(chǎn)生的離心力和力偶共同構(gòu)成飛輪微振動的兩個主要振源因素，可能激勵飛輪的模態(tài)和加劇轉(zhuǎn)子渦動，從而引發(fā)結(jié)構(gòu)共振。飛輪轉(zhuǎn)子高速旋轉(zhuǎn)時產(chǎn)生的寬頻微振動通過飛輪基座傳遞至航天器艙板，成為影響航天器姿態(tài)穩(wěn)定和指向精度的最大干擾源[5-6]。隨著用戶對航天器功能性能需求的提高，飛輪微振動對航天器平臺穩(wěn)定性的影響日益突顯出來，制約著載荷性能的進(jìn)一步提升，因此針對飛輪微振動的動力學(xué)特性研究[4]和抑制[7-8]成為慣性技術(shù)發(fā)展的重要方向。

文獻(xiàn)[9]分析了轉(zhuǎn)子振動加速度信號的頻譜，指出飛輪徑向振動頻譜中主要包括工頻及高次倍頻成分。文獻(xiàn)[10]等研究了飛輪的動力學(xué)模型，將飛輪的振源簡化為離散的諧波。文獻(xiàn)[11-12]擬合了飛輪的多階諧波。可見對飛輪微振動研究主要集中于飛輪基座上輸出的微振動頻率成分的來源定位和幅值擬合。它們是振源經(jīng)過飛輪結(jié)構(gòu)傳遞的結(jié)果，掩蓋了振源信息?？墒秋w輪振源幅值信息是飛輪產(chǎn)品微振動抑制和超靜性能提升的基礎(chǔ)，有必要對其進(jìn)行深入研究。比較簡便有效的方法是通過統(tǒng)計已有的飛輪微振動測試數(shù)據(jù)，提取飛輪振源幅值特征參數(shù)。

通過將飛輪固定在剛度較大的六分量微振動測量臺上測試飛輪振動力和力矩(如圖2所示)，忽略測試系統(tǒng)撓性對飛輪中、低頻微振動幅值的影響。測得某50 Nms飛輪的微振動力/力矩瀑布圖如圖3所示。其中，橫軸為頻率，縱軸為轉(zhuǎn)速，豎軸為振幅。

由圖3可知，飛輪微振動輸出包含了一系列放射狀倍頻振動，當(dāng)它們與轉(zhuǎn)子渦動曲線和模態(tài)曲線耦合時，可能產(chǎn)生較大結(jié)構(gòu)共振，使得振源幅值在不同的轉(zhuǎn)速和頻率點處存在不同的放大。為了從這些振動數(shù)據(jù)中提取飛輪振源特征參數(shù)，并研究其分布特性，本文首先利用魯棒回歸從25個50Nms飛輪的微振動瀑布圖中辨識出各飛輪各階次(即：微振動頻率對旋轉(zhuǎn)頻率之比)的振幅對轉(zhuǎn)速二次冪系數(shù)，將其作為該飛輪振源幅值特征參數(shù)；其次，對多個飛輪的同一階次特征參數(shù)進(jìn)行Weibull分布擬合，并計算一些近似均勻的概率點的分位數(shù)，它們既宏觀地展示了飛輪振源幅值的分布特性，也可將其作為飛輪振源幅值的分級標(biāo)準(zhǔn)。期望能為飛輪不平衡量和軸承預(yù)緊力波動量的評估和篩選提供參考。

1 基于魯棒回歸的特征參數(shù)提取

通過對飛輪微振動頻率成分進(jìn)行理論分析和測試統(tǒng)計，結(jié)果顯示飛輪的微振動階次[4]固定，其微振動特性的差異主要是振幅差異，則針對飛輪主要階次的微振動幅值統(tǒng)計具有可行性。而輸出的振幅是振源幅值經(jīng)過飛輪內(nèi)部結(jié)構(gòu)傳遞作用的結(jié)果。從傳遞角度分析，未激發(fā)結(jié)構(gòu)共振的低頻振動在結(jié)構(gòu)中的振動傳遞率近似為1；激發(fā)結(jié)構(gòu)共振的中頻段將被放大；而高頻段存在一定的衰減，這表明：

1)對于未激發(fā)結(jié)構(gòu)共振的頻率成分的振幅主要取決于飛輪不平衡量和軸承的滾動面波紋等因素的激勵作用，它反映了飛輪在該階次的激勵幅值。鑒于飛輪轉(zhuǎn)子在高速旋轉(zhuǎn)過程中產(chǎn)生的各階次微振動幅值是轉(zhuǎn)速的二次冪函數(shù)[10-13]，不妨假設(shè)第j階振幅Fj(Ω)=βjΩ2?？蓪⒃撓禂?shù)βj作為該飛輪在該階次振源幅值的一個特征參數(shù)，它較好地表征了飛輪內(nèi)部的激勵幅值隨轉(zhuǎn)速的變化規(guī)律。因此只需辨識出飛輪的主要階次的特征參數(shù)，即可比較準(zhǔn)確地反映出飛輪的微振動激勵特征。

2)低頻振動更真實反映了飛輪的微振動激勵特征，因此在進(jìn)行數(shù)據(jù)統(tǒng)計時應(yīng)賦予其更大的權(quán)重，削弱其對共振處微振動數(shù)據(jù)的依賴。鑒于普通最小二乘法(OLS)易受離群值(共振點的振動數(shù)據(jù))的杠桿作用使得對數(shù)據(jù)的擬合欠佳，擬合得到的模型的解釋性和預(yù)測性相對較差，而M估計[14]、中位數(shù)回歸和魯棒回歸[15-16]等方法都具有較好的魯棒性，對離群值不敏感。因此本文將采用魯棒回歸提取飛輪在主要頻率點處的特征參數(shù)。

對于線性歸回，不妨假設(shè)待估參量βk是n維列向量，Xk是mk×n觀測矩陣，Yk是mk維觀測列向量，Ek是mk維殘差列向量，且Yk=Xkβk+Ek，則第k次迭代的殘差加權(quán)(權(quán)矩陣為Λk)平方和：

(1)

由此可解得估計系數(shù)為：

(2)

若式(2)中Λk=I，則變?yōu)槠胀ㄗ钚《嘶貧w：

(3)

則觀測向量Yk的OLS估計值為：

(4)

(5)

依據(jù)殘差可計算出殘差中位數(shù)Mk為：

(6)

式中：median為中位數(shù)函數(shù)。

誤差項的標(biāo)準(zhǔn)離差估計sk：

(7)

則其標(biāo)準(zhǔn)化殘差rk為：

(8)

其中，εk為前一次回歸的殘差；tc為調(diào)節(jié)常數(shù)。若將加權(quán)最小二乘的權(quán)重取為殘差的函數(shù)(對殘差較大的數(shù)據(jù)賦予較小的權(quán)重，反之亦然)，其變?yōu)轸敯艋貧w。常用的權(quán)函數(shù)有Bisquare核函數(shù)、Cauchy核函數(shù)和Welsch核函數(shù)等，其中，Bisquare核函數(shù)：

(9)

(10)

計算流程如圖4所示。由此可回歸得到飛輪各階次振源幅值的二次冪系數(shù)βj(j∈{1,2,…,∞}),將其作為飛輪的振源特征參數(shù)。當(dāng)飛輪結(jié)構(gòu)不變時,系數(shù)βj也將不變，它是飛輪的不變特征量。

2 基于Weibull分布的特征參數(shù)分級

50 Nms飛輪產(chǎn)品的生產(chǎn)線和工藝過程固化，產(chǎn)品的研制質(zhì)量嚴(yán)格受控，其微振動特性參數(shù)相對穩(wěn)定，所積累的多個產(chǎn)品的微振動數(shù)據(jù)具有一定的統(tǒng)計規(guī)律性，故可通過統(tǒng)計的方法研究飛輪的振源特征參數(shù)分布特性。但微振動特性受到零件參數(shù)、裝配間隙、軸承預(yù)緊力等眾多因素的影響，因此又很難從理論上證明飛輪的微振動特征參數(shù)一定服從某種概率分布，故需進(jìn)行試探性擬合。鑒于飛輪特征參數(shù)都大于零，根據(jù)直方圖可見，其分布形式與廣義Gamma分布或者Weibull分布[17]近似。該分布參數(shù)的取值可通過擬合得到，可將其作為該飛輪的分布特征參數(shù)。廣義Gamma分布的概率密度函數(shù)：

(11)

式中:x>0，ρ>0為尺度參數(shù)；α>0為功率參數(shù)；λ>0為形狀參數(shù)。Г(λ)為Gamma函數(shù)。當(dāng)λ=1時，即可得到Weibull分布的概率密度和分布函數(shù)：

(12)

(13)

利用該Weibull分布來擬合各階次微振動幅值特征參數(shù)，求出分布參數(shù)。則其逆函數(shù)為：

(14)

為了研究飛輪振源特征參數(shù)的分布特性，不妨以近似均勻的概率點p∈{0.2,0.4,0.6,0.8,0.95}作為其微振動分級概率，飛輪等級與分位數(shù)的對應(yīng)關(guān)系如表1所示,再用式(14)求出表1的分位數(shù)。

表1 微振動特性分級表Table 1 Rating table of micro vibration characteristics

3 數(shù)據(jù)統(tǒng)計結(jié)果

飛輪主要頻率階次統(tǒng)計結(jié)果列于表2第2列。針對主要階次在微振動瀑布圖(緩慢升速時的振動頻譜曲線隨轉(zhuǎn)速變化形成的三維圖)中搜索振動幅值隨轉(zhuǎn)速變化的曲線，再利用圖4求得該飛輪該階次的二次特征參數(shù)βj。如某飛輪的工頻和3倍頻在升速過程中的振幅曲線及其擬合曲線如圖5所示。

圖5中曲線上振動較大的尖峰主要是因為該階次振源激勵被飛輪的結(jié)構(gòu)共振放大所致。對未激發(fā)結(jié)構(gòu)共振的工頻振動的幅-頻曲線具有較好的擬合精度，多次擬合曲線基本重合；對于存在結(jié)構(gòu)共振的幅-頻曲線，需采用圖4的計算流程逐步消除結(jié)構(gòu)共振導(dǎo)致的離群值，擬合曲線終將收斂于該階次的激勵幅值。該擬合曲線在一定程度上直觀反映了該飛輪的微振動激勵幅值隨轉(zhuǎn)速變化的關(guān)系。

利用圖4的計算流程逐一提取飛輪各主要階次的微振動激勵幅值特征參數(shù)，進(jìn)而可得到飛輪在全轉(zhuǎn)速范圍內(nèi)的激勵幅值-頻率-轉(zhuǎn)速三維圖(某飛輪的激勵三維圖如圖6所示)。與振動輸出(圖3所示)對比可見，回歸得到的激勵三維圖消除了結(jié)構(gòu)共振影響，直觀反映了飛輪微振動激勵幅值對頻率和轉(zhuǎn)速的關(guān)系。結(jié)果表明該飛輪除了工頻(不平衡量導(dǎo)致的離心力/力偶)幅值較大以外，其余由軸承激發(fā)的預(yù)緊力波動量的幅值較小?？梢姡w輪中、高頻振動主要取決于這些微小振動被結(jié)構(gòu)共振放大的倍數(shù)。

從25個50 Nms飛輪的微振動瀑布圖數(shù)據(jù)中逐一提取各飛輪各主要階次的振源激勵特征參數(shù)，如表2所示。對表2中各階次的所有特征參數(shù)按照Weibull分布進(jìn)行擬合，從而得到各主要階次的分布參數(shù)，再根據(jù)式(14)計算表1要求的分位數(shù)?？傻玫礁髦饕A次微振動激勵幅值的分級圖如圖7所示，圖中將飛輪的振動激勵幅值分為六個等級，其中Ⅵ級為超出Ⅴ級以上的廣泛區(qū)域。

從表2及圖6和圖7可以看出，飛輪的不平衡特征參數(shù)較大，且分布范圍較寬，主要位于10-8量級；而飛輪軸承預(yù)緊力波動量的主要階次的特征參數(shù)主要位于10-9量級。除{0.6,1,1.78,2.23,3.43,5.59}倍頻的振動較大以外，其它倍頻振動相對較小，但由于飛輪的結(jié)構(gòu)阻尼也較小，這些微小的倍頻振動也可能會激發(fā)結(jié)構(gòu)共振，從而輸出較大的振動。因此抑制飛輪的振動可主要從如下三個方面著手：1)減小飛輪的剩余不平衡量；2)增加飛輪的結(jié)構(gòu)阻尼；3)研究飛輪的共振機(jī)理，優(yōu)化模態(tài)曲線和渦動曲線分布，使之避開振幅較大的激勵。

表2 25個飛輪在主要階次的特征參數(shù)統(tǒng)計表Table 2 Characteristic parameters statistics at some main orders of 25 flywheels

4 結(jié) 論

為了研究飛輪振源幅值的分布特性，本文利用魯棒回歸的方法提取了與振源幅值直接相關(guān)的轉(zhuǎn)速二次冪系數(shù)，將其作為飛輪的振源特征參數(shù)；并運(yùn)用統(tǒng)計的方法對同階次的特征參數(shù)進(jìn)行擬合，計算出特定概率點的分位數(shù)。將該方法成功應(yīng)用于對50 Nms飛輪的微振動數(shù)據(jù)統(tǒng)計分析中，從微振動瀑布圖中提取了飛輪振源幅值特征參數(shù)，所計算出的分位數(shù)直觀反映了該型飛輪的特征參數(shù)的分布情況。統(tǒng)計結(jié)果亦可作為該型飛輪振源特征參數(shù)的分級標(biāo)準(zhǔn)，為超靜飛輪的篩選提供初步的判定準(zhǔn)則，為飛輪改進(jìn)提供參考依據(jù)，為飛輪在某種試驗前后的結(jié)構(gòu)一致性比對提供一種簡單有效的方法。

[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] 王虹,卿濤,周剛,等.預(yù)載對飛輪用軸承組件振動特性的影響研究[J].空間控制技術(shù)與應(yīng)用, 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的高速轉(zhuǎn)子的徑向振動特性研究[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] 劉天雄,范本堯,楊慧.衛(wèi)星飛輪擾振控制技術(shù)研究[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] 龐世偉,楊雷,曲廣吉.高精度航天器微振動建模與評估技術(shù)最近進(jìn)展[J].強(qiáng)度與環(huán)境,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] 劉天雄,林益明,王明宇,等.航天器振動控制技術(shù)進(jìn)展[J]. 宇航學(xué)報, 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]. 宇航學(xué)報, 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].空間控制技術(shù)與應(yīng)用, 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] 趙蘭敏,郄志紅,吳鑫淼,等.魯棒回歸在水工建筑物安全監(jiān)控模型建立中的應(yīng)用[J].中國農(nóng)村水利水電,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初值的穩(wěn)健回歸的諧波發(fā)射水平評估方法[J].電力系統(tǒng)保護(hù)與控制, 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分布參數(shù)的收縮估計[J].宇航學(xué)報, 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.]