連續(xù)轉(zhuǎn)速狀態(tài)下衛(wèi)星反作用飛輪微振動參數(shù)識別

楊 林 ，王巖松 *，魏 磊 ，胡自強

（1. 山東大學前沿交叉科學青島研究院，山東青島266200；2. 山東大學 空間科學研究院，山東 威海 264209）

1 引 言

隨著航天技術(shù)的快速發(fā)展，光學遙感衛(wèi)星正呈現(xiàn)出更高指向精度和更高分辨率的特點。近些年國際上的商業(yè)遙感衛(wèi)星已經(jīng)達到優(yōu)于0.3 m的分辨率［1］。對于高分辨率遙感衛(wèi)星來說，光學載荷對于衛(wèi)星上的微小擾動十分敏感，微振動對成像質(zhì)量的影響已經(jīng)變得不可忽略。在反作用飛輪微振動作用下，光學遙感衛(wèi)星在成像時會出現(xiàn)圖像扭曲和模糊等問題［2］。

為降低衛(wèi)星微振動的影響，主要采取的手段包括：降低干擾源的影響、優(yōu)化星體傳力路徑以及隔振等方法［3-5］。反作用飛輪是航天器的姿態(tài)控制系統(tǒng)的執(zhí)行機構(gòu)之一，其性能影響著航天器姿態(tài)控制的精度。由于反作用飛輪在工作時受到轉(zhuǎn)子動靜不平衡、軸承缺陷、飛輪結(jié)構(gòu)模態(tài)放大等因素的影響［6］，會激發(fā)出一系列復雜的諧波擾動，所以反作用飛輪是微振動的主要干擾源之一。對反作用飛輪力學機理以及其運轉(zhuǎn)過程中的干擾特性的研究，對航天器微振動技術(shù)以及其抑制方法的研究具有重要的意義。對反作用飛輪的研究方法主要包括飛輪理論模型的建模［7］，飛輪干擾特性的經(jīng)驗模型的建立以及飛輪干擾特性及結(jié)構(gòu)模態(tài)參數(shù)識別等［8-11］。通過飛輪參數(shù)識別方法可以直觀地反應出飛輪干擾特性和結(jié)構(gòu)模態(tài)的分布特點，并且能夠反應出理論建模不能反應的真實特性，為后續(xù)飛輪的隔振方法研究、整星的微振動分析和試驗提供基礎(chǔ)。

當前應用范圍較廣的是恒定轉(zhuǎn)速的飛輪參數(shù)識別方法，通常是利用一系列固定轉(zhuǎn)速的干擾力/力矩數(shù)據(jù)，分析飛輪的干擾特性及模態(tài)參數(shù)［12］。此類方法僅能分析恒定轉(zhuǎn)速下飛輪的模態(tài)參數(shù)和干擾力特性，為了得到飛輪在一個轉(zhuǎn)速區(qū)間的干擾特性，需要通過大量的試驗來獲得不同轉(zhuǎn)速情況下試驗數(shù)據(jù)［13］，大大增加了工作量；并且，反作用飛輪在實際工作過程中，尤其在衛(wèi)星姿態(tài)運動過程中，其轉(zhuǎn)速并不是完全恒定的［14］，飛輪的干擾力為一段非平穩(wěn)信號，此時利用傳統(tǒng)的傅里葉變換方法得到的模態(tài)特性和干擾力特性將因轉(zhuǎn)速的變化而變得不穩(wěn)定。所以常規(guī)的恒轉(zhuǎn)速下飛輪微振動參數(shù)識別方法是不能準確反應在軌工作狀態(tài)的。

本文將建立一種非恒定轉(zhuǎn)速運動狀態(tài)下反作用飛輪的參數(shù)識別方法，可以利用非平穩(wěn)信號識別出飛輪非恒定轉(zhuǎn)速下的干擾力倍頻特性以及模態(tài)參數(shù)。首先利用時頻分析，將試驗得到的飛輪非恒定轉(zhuǎn)速下的干擾力非平穩(wěn)信號，轉(zhuǎn)化為時頻分布函數(shù)，并通過轉(zhuǎn)速與時間的對應關(guān)系，將時頻分布函數(shù)轉(zhuǎn)化為轉(zhuǎn)頻分布函數(shù)；之后建立飛輪的倍頻特性識別方法，識別出飛輪的干擾力分布；最后建立非恒定轉(zhuǎn)速飛輪的參數(shù)化的結(jié)構(gòu)模態(tài)參數(shù)識別方法，識別出飛輪的結(jié)構(gòu)模態(tài)。本文將通過干擾力測試驗證所提出的方法，首先進行一次非恒定轉(zhuǎn)速的干擾力測試得到飛輪干擾力試驗數(shù)據(jù)，利用以上方法識別出該飛輪的干擾力分布特性和模態(tài)參數(shù)。通過結(jié)合倍頻特性和模態(tài)參數(shù)，可以得到飛輪的轉(zhuǎn)速-頻率峰值點，預測飛輪干擾力最大值的分布情況。該方法將在保證一定精度的前提下，大大提升飛輪參數(shù)識別的效率。

2 反作用飛輪倍頻特性識別

2.1 動力學模型

反作用飛輪的動力學模型可以表示為

通過求解特征方程，可以得到反作用飛輪的無阻尼固有頻率，可以表示為［13］：

其中：ωr和ωa分別對應徑向平動和軸向平動，ωψ，?對應徑向搖擺，Ω指飛輪轉(zhuǎn)速，ξ指極慣性矩與徑向慣性矩之比。從式（2）中可以發(fā)現(xiàn)，飛輪的搖擺模態(tài)是與飛輪轉(zhuǎn)子速度相關(guān)的。

2.2 時頻分析原理

對于固定轉(zhuǎn)速的飛輪，其干擾力/力矩時間響應信號是一組平穩(wěn)信號，通過Fourier 變換即可得到當前轉(zhuǎn)速下飛輪各個方向干擾力的頻譜。但當飛輪的轉(zhuǎn)速不斷變化時，飛輪干擾力信號的頻譜特性是隨時間變化的。因此傳統(tǒng)的Fourier變換不再適應變轉(zhuǎn)速飛輪干擾力/力矩頻域特性分析。

時頻分析方法是一種時變非平穩(wěn)信號處理方法，能夠?qū)⒎瞧椒€(wěn)信號分解為由時間和頻率二元變量決定的時頻分布函數(shù)，以描述信號在不同時間和頻率下的幅值分布密度［15］。時頻分析方法可以直觀的體現(xiàn)信號的頻率變化特性，廣泛應用于故障診斷和模態(tài)參數(shù)識別等領(lǐng)域。目前應用最為廣泛的時頻分析方法主要包括短時Fourier 變換（STFT）、小波變換（WT）、Wigner-Ville分布（WVD）等。以下將采用短時Fourier 變換作為飛輪響應信號的時頻分析方法。

短時Fourier 變換的頻譜可以表示為：

其中：x(τ)為被分析的時間信號，h(t)為窗函數(shù)，ht，ω(τ) =h(τ-t)e-j2πfτ為窗函數(shù)對應的基函數(shù)。

之后利用轉(zhuǎn)速隨著時間變化的函數(shù)，可以將時頻分布譜的時間軸轉(zhuǎn)化為轉(zhuǎn)速軸，得到短時Fourier 變換在轉(zhuǎn)速- 頻率平面上的分布STFT(r，f)，功 率 譜 密 度 函 數(shù) 可 以 表 示 為G(r，f) =|STFT(r，f)|2。

2.3 倍頻特性識別

飛輪的干擾力由一系列復雜的諧波和噪聲構(gòu)成。這些諧波包括基礎(chǔ)諧波、超諧波和次諧波，其中基礎(chǔ)諧波即為一階諧波，而超諧波和次諧波對應的頻率均可以表示為一階諧波頻率的倍數(shù)，即倍頻。本節(jié)將通過建立線性最小二乘方法，識別出飛輪各個方向干擾力或力矩的主要倍頻。

定義飛輪倍頻特性的線性最小二乘估計的費用函數(shù)為：

其中：系數(shù)p為待估計的倍頻系數(shù)，ri為第i個轉(zhuǎn)頻采樣點的轉(zhuǎn)速，fi為第i個轉(zhuǎn)頻采樣點的頻率值，Nr和Nf分別為轉(zhuǎn)速點數(shù)和頻率點數(shù)，誤差函數(shù)可以表示為：

其中：wi為第i個時頻采樣點的權(quán)值，引入一個如圖1 所示的扇形窗函數(shù)，窗函數(shù)可以表示為：

其中D定義了一個域，在域中p1≤r/f＜p2。

將第2.2 節(jié)得到的時頻分布函數(shù)G(r，f)寫為向量形式：

圖1 扇形窗Fig.1 Sector window

并令權(quán)函數(shù)wi=hsiGi，則誤差函數(shù)可以表示為：

其中：符號.*表示兩向量元素相乘，待估參數(shù)p的最小二乘解可以表示為

通過改變扇形窗函數(shù)的分布，即可識別出不同扇形窗內(nèi)的倍頻特性。

3 反作用飛輪模態(tài)參數(shù)識別

3.1 參數(shù)化的矩陣分式模型建立

本節(jié)將建立一種與轉(zhuǎn)速相關(guān)的參數(shù)化的飛輪結(jié)構(gòu)模態(tài)參數(shù)識別方法，首先建立飛輪的參數(shù)化矩陣分式模型［16］，利用加權(quán)非線性最小二乘方法估計出模型的待估參數(shù)，最后將矩陣分式模型參數(shù)轉(zhuǎn)化為模態(tài)參數(shù)。

飛輪的動力學模型可以寫為如下形式的參數(shù)化的右矩陣分式模型：

其中：Bk為分子系數(shù)矩陣，k表示第k個輸出信號，A為分母系數(shù)矩陣，r表示轉(zhuǎn)速變量，f表示頻率變量。對于單參考情況，可以簡化為公分母模型：

將分母系數(shù)和分子系數(shù)定義為如式（12）所示的多項式形式：

其中：φi，j(r，f)表示與頻率和轉(zhuǎn)速相關(guān)的基函數(shù)，采用正交多項式與z域多項式混合形式φi，j(r，ω) =pi(r)zj，bk，i，j和ai，j為 多 項 式 系 數(shù) ，nr和nf分別為轉(zhuǎn)速多項式和頻率多項式的階數(shù)［17］。

將多項式系數(shù)寫成向量形式如下：

則右矩陣分式模型的所有待估參數(shù)寫成以下形式：

3.2 加權(quán)非線性最小二乘法處理

本小節(jié)利用加權(quán)非線性最小二乘估計方法辨識出參數(shù)化模型的待估參數(shù)，將模型參數(shù)轉(zhuǎn)化為模態(tài)參數(shù)。加權(quán)非線性最小二乘估計的費用函數(shù)可以表示為：

式中，符號H 表示矩陣的共軛轉(zhuǎn)置，第k個輸出信號的誤差函數(shù)εk(rl，fm，θ)可以表示為：

加權(quán)非線性最小二乘方法的待估參數(shù)θ需要通過迭代的方式計算，采用Gauss-Newton 算法，待估參數(shù)θ的增量 Δθ可由如式（21）的線性方程進行求解：

其中：J(θn)表示第n個迭代步的 Jacobi矩陣，Jacobi矩陣J的含義是誤差函數(shù)對待估參數(shù)θ的偏導：

Jacobi 矩陣中對應的分塊矩陣如式（23）所示：

其中?表示kronecker 積，

則式（21）可以表示為：

根據(jù)式（26）的第No+1 行，可得：

結(jié)合式（27）和式（28）消去 Δβk，可得到縮減正則方程如式（29）：

定義矩陣：

則通過求解式（31）的方程即可得到每一次迭代的 Δα和 Δβk：

將每步得到的 Δα和 Δβk合寫成 Δθ，并對待估參數(shù)θ進行更新，即：

迭代終止條件為：

其中：ζ為相對誤差的閾值，kiter為迭代次數(shù)，kmax為給定的最大迭代次數(shù)。

待估參數(shù)的初始值θ0可由線性最小二乘估計方法所估計出的參數(shù)給定［18］。

3.3 模態(tài)參數(shù)計算

模態(tài)參數(shù)可以通過求解如式（34）廣義伴隨矩陣的特征值獲得：

則模態(tài)頻率fr和阻尼比ξr可以分別表示如式（38）所示：

其中：λr為廣義伴隨矩陣Ac(rj)的特征值，即λr=eig(Ac(rl))。

4 試驗與分析

本節(jié)將通過連續(xù)轉(zhuǎn)速狀態(tài)下的反作用飛輪干擾力/力矩測試，來驗證本文提出的方法。同時增加一組反作用飛輪固定轉(zhuǎn)速下的干擾力/力矩測試，通過傳統(tǒng)的傅里葉變換方法進行分析，并與本文所提出方法的結(jié)果進行對比。

根據(jù)第2 節(jié)和第3 節(jié)內(nèi)容，編寫分析程序，程序偽代碼如圖2 所示。

4.1 試驗工況

如圖3 所示，試驗系統(tǒng)由飛輪、kistler 測力臺、試驗工裝、信號采集系統(tǒng)、控制器、電源和計算機組成。將飛輪通過試驗工裝連接在測力臺上，飛輪方向如圖4 所示定義。設(shè)定飛輪力矩為0.5 N·m，測試飛輪轉(zhuǎn)速從0 轉(zhuǎn)勻加速轉(zhuǎn)到2 000轉(zhuǎn)的過程，在測得干擾力信號的同時，記錄飛輪轉(zhuǎn)速隨時間變化的曲線。測試的現(xiàn)場圖如圖5所示。

試驗得到的飛輪的轉(zhuǎn)速時間歷程曲線如圖6所示，飛輪各個方向干擾力/力矩的時間歷程曲線如圖7 所示。可以發(fā)現(xiàn)，連續(xù)轉(zhuǎn)速狀態(tài)下飛輪的干擾力信號是一條幅值呈上升趨勢的非穩(wěn)態(tài)信號。此時利用傳統(tǒng)的傅里葉變換方法進行分析時，干擾力在頻域上的頻率和幅值都將是不穩(wěn)定的，因此需要將連續(xù)轉(zhuǎn)速的干擾力信號轉(zhuǎn)化的時頻域上進行分析。

4.2 干擾力/力矩時頻分析

利用第2.2 節(jié)的方法將干擾力信號轉(zhuǎn)化到時頻域，得到干擾力和力矩的轉(zhuǎn)頻譜如圖8 所示。圖中可以發(fā)現(xiàn)明顯的飛輪結(jié)構(gòu)模態(tài)和倍頻特性，可以通過本文建立的方法進行識別。

圖2 分析程序偽代碼Fig.2 Pseudo-code of analysis program

圖3 試驗系統(tǒng)示意圖Fig.3 Schematic diagram of experimental system

圖4 飛輪安裝方向坐標定義Fig.4 Flywheel installation direction definition

圖5 測試現(xiàn)場圖Fig.5 Photo of experimental site

圖6 轉(zhuǎn)速時間歷程曲線Fig.6 Time history of rotating speed

圖7 干擾力和力矩信號的時間歷程曲線Fig.7 Time responses of the disturbing forces/moments signals

4.3 倍頻特性識別

利用第2.3 節(jié)建立的方法對飛輪倍頻特性進行識別，識別出的各方向干擾力或力矩的倍頻特性進行統(tǒng)計如表1 和圖9 所示。

4.4 飛輪模態(tài)參數(shù)識別

利用第3 節(jié)建立的方法辨識飛輪的模態(tài)參數(shù)，設(shè)置轉(zhuǎn)速多項式參數(shù)nr= 1，頻率多項式參數(shù)，最大迭代次數(shù)kmax= 10，辨識出的轉(zhuǎn)頻譜如圖10 所示。辨識的隨轉(zhuǎn)速變化的模態(tài)頻率和阻尼比如圖11 所示，圖中，紅色矩形符號表示飛輪的軸向平動模態(tài)，綠色圓形符號表示飛輪的搖擺模態(tài)，藍色菱形符號表示飛輪的徑向平動模態(tài)（彩圖見期刊電子版）。所辨識出的飛輪模態(tài)可以為飛輪的仿真分析、模型校核以及飛輪系統(tǒng)隔振設(shè)計等提供參考。

圖8 干擾力和力矩的轉(zhuǎn)頻譜Fig.8 Speed-frequency spectra of forces and moments

圖9 倍頻辨識結(jié)果（各直線表示各倍頻辨識結(jié)果）Fig.9 Identification results of harmonic characteristics（Lines denote the identification results）

表1 飛輪各方向干擾力/力矩倍頻統(tǒng)計Tab.1 Summary of main harmonic factors

圖10 參數(shù)化估計的轉(zhuǎn)頻譜Fig.10 Parametric estimated speed-frequency spectra

圖11 模態(tài)頻率和阻尼比Fig. 11 Modal frequency and damping ratio

4.5 固定轉(zhuǎn)速試驗結(jié)果

在進行變轉(zhuǎn)速飛輪干擾力測試的同時，進行固定轉(zhuǎn)速的干擾力測試作為對比試驗。從0 轉(zhuǎn)到2 000 轉(zhuǎn)每隔100 轉(zhuǎn)進行一組測試，共進行21 組試驗。通過傅里葉變換，得到飛輪的干擾力/力矩隨頻率和轉(zhuǎn)速變化的瀑布圖如圖12 所示。

通過結(jié)合飛輪模態(tài)頻率和倍頻特性，可以預測干擾力或力矩的峰值點位置。如圖13 所示（彩圖見期刊電子版），將飛輪軸向模態(tài)頻率和倍頻曲線顯示在同一張圖上，兩曲線交點即為干擾力轉(zhuǎn)頻峰值點。通過表2 給出了連續(xù)轉(zhuǎn)速試驗和固定轉(zhuǎn)速試驗的軸向干擾力的轉(zhuǎn)頻峰值點及其相對誤差，可以發(fā)現(xiàn)兩種方式得到的峰值點誤差最大為4.3%，平均誤差經(jīng)計算小于1%。并且固定轉(zhuǎn)速試驗得到的峰值點只能出現(xiàn)在試驗測試設(shè)定的轉(zhuǎn)速上，容易忽略最大峰值點所出現(xiàn)的轉(zhuǎn)速。而連續(xù)轉(zhuǎn)速方法可以得到干擾力（或力矩）的平滑變化趨勢。此外，連續(xù)轉(zhuǎn)速試驗的試驗量僅為固定轉(zhuǎn)速試驗的4.8%。

圖12 固定轉(zhuǎn)速試驗瀑布圖Fig.12 Waterfall plots of the fixed-speed tests

圖13 軸向干擾力峰值分布（紅色實線表示模態(tài)頻率，黑色虛線表示倍頻曲線，藍色點表示參考試驗的峰值點）Fig.13 Distribution map of the peak points of Fz（red line denotes modal frequencies，black dash lines are the harmonic frequencies，blue points denote the peak points）

表2 軸向干擾力的轉(zhuǎn)頻峰值點Tab.2 Speed-frequency points of peaks of Fz

5 結(jié) 論

本文建立了一種連續(xù)轉(zhuǎn)速下反作用飛輪的參數(shù)識別方法，該方法利用連續(xù)轉(zhuǎn)速干擾力數(shù)據(jù)，識別出飛輪的模態(tài)參數(shù)和倍頻特性，所識別的結(jié)果相對于固定轉(zhuǎn)速參考試驗的結(jié)果基本一致，但是卻獲得了固定轉(zhuǎn)速測試所不能得到的干擾力（力矩）平滑變化趨勢數(shù)據(jù)，對飛輪擾振特性有更準確地了解。通過結(jié)合倍頻特性和模態(tài)頻率特性，可以計算出飛輪干擾力和力矩的轉(zhuǎn)頻峰值點。本文提出的方法提高了飛輪參數(shù)辨識的效率，揭示了飛輪干擾力/力矩隨轉(zhuǎn)速變化的關(guān)系，可以用于反作用飛輪特性的獲取，識別的模態(tài)頻率、阻尼比、倍頻以及轉(zhuǎn)頻峰值點等特性，可以為后續(xù)飛輪模型校核、仿真分析以及整星微振動分析與試驗等微振動相關(guān)技術(shù)的研究提供參考。