基于SGHT 的航天器微振動信號處理方法?

陳繼喆， 章海鷹

（1.中國科學院國家天文臺南京天文光學技術研究所 南京，210042）

（2.中國科學院天文光學技術重點實驗室（南京天文光學技術研究所） 南京，210042）

（3.中國科學院大學 北京，100049）

（4.中國科學院大學天文與空間科學學院 北京，100049）

引 言

航天器飛輪微振動是影響空間光學載荷成像質量的主要因素［1］。為了評估微振動的影響，在研制階段需要開展地面微振動試驗［2］。常見的試驗方法是多點布置壓電振動傳感器，采集飛輪工作時不同位置的加速度信號。這種方法測量難度低，且可以對航天器內部進行測量，但得到的航天器微振動加速度信號需要進一步處理［3］。

通常情況下，希望經過處理后的微振動信號提供兩方面的信息：①頻域上分析微振動能量集中的特定頻段，指導載荷結構優化以避免與飛輪微振動發生共振；②時域上通過對微振動加速度信號二次積分，得到微振動線位移、角位移估計，判斷平臺擾動對光學成像的具體影響。因此，研究如何處理航天器微振動信號是開展地面微振動試驗的必要保障。

為得到頻域信息，考慮微振動信號的非平穩特性，可以采用希爾伯特-黃轉換（Hilbert-Huang transform，簡稱HHT）進行處理，即將原始信號經過經驗模態分解（empirical mode decomposition，簡稱EMD）處理后，對各信號分量通過希爾伯特轉換（Hilbert transform，簡稱HT）求解瞬時頻率，繪制時間-頻率-瞬時能量瀑布圖，得到包含頻域能量分布特征的邊際譜。這種信號處理方法簡單，但缺點是EMD分解過程中容易產生欠包絡、過包絡和頻率混淆等問題，導致得到的頻率成分可能在實際信號中并不存在，從而影響對微振動特征的識別。

為得到時域信息，常見的方法包括頻域積分法和時域積分法。頻域積分法通過快速傅里葉變換（fast Fourier transform，簡稱FFT）在頻域實現加速度與位移的相互轉換［4］；時域積分法通過例如辛普森公式的數值積分方法，直接對時域加速度信號進行積分。這兩種方法都可以得到時域位移數據，但都存在低頻段微小噪聲和誤差引發二次積分后位移出現巨大趨勢偏差的問題。為了消減趨勢偏差影響，目前有多種常用方法可以對加速度信號和積分位移進行處理。

從加速度信號入手，考慮到低頻段噪聲誤差是引起位移趨勢偏差的主要原因，高通濾波算法［5］和低頻衰減濾波算法［6］選擇舍去部分低頻段加速度信號特征，對趨勢偏差有一定改善，但需要根據工程經驗設置低頻參數，且沒有消除通帶內信號的噪聲，噪聲魯棒性較差。小波積分方法盡管曾用于處理量子科學實驗衛星的微振動加速度信號［3］，但同樣需要選擇恰當的窗函數構建衰減曲線［7］，也存在低頻衰減算法的眾多缺點。

從積分位移入手，位移的長周期趨勢項是趨勢偏差的主要構成。多項式去趨勢方法通過對積分位移進行多項式擬合，將擬合的多項式剔除以去除趨勢項。該方法主要針對因為初值不確定產生的趨勢項誤差，對其他原因產生的趨勢項誤差擬合效果較差。EMD 方法分解得到的最長周期殘余信號通?？梢源碓盘柕内厔莼蚓?，因此也可以用于信號的去趨勢處理［8］，但實際工程中受限于分解能力，效果并不理想。

針對現有常用航天器微振動信號處理方法的不足，筆者提出SGHT 方法，自適應地完成頻域特征獲取及時域信號處理。首先，對輸入的時域信號進行SGMD［9-10］，得到各信號分量，具有良好的噪聲魯棒性和分解性能；其次，引入歸一化互信息作為評價指標，同時參考迭代式SGMD 方法［11］，迭代式的輸出單信號分量歸類重組的SGC，進一步增強自適應能力、噪聲魯棒性和分解性能；最后，應用希爾伯特譜分析，得到各信號分量的瞬時頻率，并以頻域能量分布的形式輸出結果。

在頻域特征獲取方面，SGHT 方法將微振動加速度信號作為原始輸入信號，可以在有噪聲干擾的情況下，分解、重構信號并輸出各SGC 分量及對應的頻譜特征，對于常用HHT 方法難以處理的含噪非平穩微振動信號，SGHT 方法具有顯著優勢。

在時域信息處理方面，首先，使用SGHT 方法處理微振動加速度信號；其次，對處理后的微振動加速度信號二次積分，得到微振動位移估計；然后，將得到的位移估計作為輸入信號再次進行SGHT 處理，得到各SGC 分量及對應的頻譜特征；最后，將包含趨勢項的SGC 分量去除，得到更準確的二次積分位移結果。

1 基本原理

本研究辛幾何-希爾伯特的轉換流程見圖1。

圖1 辛幾何-希爾伯特的轉換流程Fig.1 The process of SGHT

1.1 相空間重構

設輸入的原始信號時間序列為x=x1，x2，…，xn，其中n為數據長度，重構序列可以得到軌跡矩陣X為

其中：d為嵌入維數；τ為延遲時間，可取1；m=n-(d-1)τ。

嵌入維數用功率譜密度法確定［12］，首先計算信號時間序列x的功率譜密度（power spectral density，簡稱PSD），然后估計PSD 中最大峰值對應的頻率fmax，嵌入維數即確定為，其中Fs為 輸入信號的采樣頻率。為縮小計算規模，嵌入維數最大限定為200。

1.2 辛幾何矩陣變換

令協方差對稱矩陣A=XTX，構造Hamilton 矩陣M為

令矩陣N=M2，則N也為Hamilton 矩陣。構造正交辛矩陣Q

其中：Q為正交辛矩陣，在變換后依然保留Hamilton矩陣的結構形式；B為上三角矩陣，即矩陣元素bij=0(i>j+1)。

矩陣B可通過對矩陣N施密特正交化得到，其特 征 值 為λ1，λ2，…，λd。根 據Hamilton 矩 陣 的 性 質可知，矩陣A的特征值σi=(i=1，2，…，d)，設Vi(i=1，2，…，d)為矩陣A對應于特征值σi的特征向量。

計算轉換系數矩陣Si=ViTXT，得到重構矩陣Zi=ViSi，Zi(i=1，2，…，d)即為初始單分量矩陣。重構軌跡矩陣Z由d個組分組成，可寫為Z=Z1+Z2+…+Zd。

1.3 對角平均化

實際得到的重構軌跡矩陣Zi為m×d矩陣，為了得到長度為n的時間序列，需要對其進行對角平均化。對于任一初始單分量矩陣Zi，定義矩陣的各元 素 為zij，其 中：1 ≤i≤d；1 ≤j≤m。 令d*=min (m，d)，m*=max (m，d)，n=m+(d-1)τ。當m

由式（4）計算得到的y1，y2，…，yn即為對Zi處理得到的一組長度為n的一維時間序列Yi。對各重構軌跡矩陣Zi依次處理，可得到原信號分解的d組長度為n的分量，即原信號x=Y=Y1+Y2+…+Yd。

1.4 單分量重構

矩陣分解得到d個單分量信號后，并非每個單分量都是獨立分量，多個不同分量可能由相同的原因引起擾振，具備相似的特征。因此，本方法將歸一化 互 信 息（normalized mutual information， 簡 稱NMI）作為評價指標，將NMI 較高的分量歸類重組。同時，測得的原始信號可能還包括測量誤差和其他因素引起的噪聲，這些分量不宜歸類重組，需要設置合適的終止條件。

計算Y1與其他分量的歸一化互信息，即

其 中：EB(A)=-∑a P(a)logP(a)，為 邊 界 熵；EJ(X，Y)=-∑x∑y P(x，y)logP(x，y)，為聯合熵。

若存在其他分量與Y1的歸一化互信息高于設定的閾值，則所有符合閾值條件的分量與Y1歸類重組為SGC 分量；若不存在此條件分量，則Y1單獨組成SGC 分量。

在原始信號中去除得到的SGC 分量，得到余量信號gh+1(t)為

計算所得殘余項和原信號之間的歸一化平均絕對 誤 差（normalized mean absolute error， 簡 稱NMAE），即

其中：gh+1(i)，x(i)分別為gh+1(t)，x(t)時程的位移采樣值；Npoint為信號的時間采樣點數。

NMAE 作為歸一化指標，其最大值為殘余項信號且是原信號自身時，則NMAE 為1；其最小值為不存在殘余項信號時，則NMAE 為0。因此，根據此范圍可以將NMAE 的默認閾值確定為0.1。 當NMAE 小于給定閾值即代表原始信號已經完成分解；否則得到的余量信號將作為新一輪分解過程的輸入信號，并分解得到新的SGC 分量。

1.5 邊際譜計算

對處理后的SGC 分量進行希爾伯特變換，即

構造解析函數

其中：ai(t)為包絡信號；θi(t)為瞬時相位；ωi(t)為瞬時角頻率。

該信號的希爾伯特譜可表示為

希爾伯特邊際譜表示為

希爾伯特邊際譜幅值是每一個頻率點幅值在時間全局上的積分［13］，從統計學角度看，幅值越大表示具有該頻率的特征分量在整個信號持續時間內某一時刻出現的可能性較高，因此能夠精確反映某一頻率是否真實存在［14］，其幅值是適合定性分析的無量綱量。最終，輸出所有SGC 分量及其對應的希爾伯特邊際譜。

2 仿真驗證

2.1 微振動模型

為了模擬實際信號在頻譜中呈現“樅”狀分布的特性，在現有航天器飛輪微振動模型［15］的基礎上，假定飛輪轉動過程中轉速和幅值變化的瞬態效應不可忽略，但其在微小時間段（1 個諧波周期）內仍保持常數。因此，在原線性穩態模型中引入變化因子，微振動模型可表示為

其中：a(t)為加速度；n為模型總階數；Ci為第i階諧波的振幅系數；frwa為飛輪轉速頻率；αij和βij分別為控制幅值和周期的隨機變化因子；j為第j個諧波周期；hi為諧波級次；t為時間；φi為［0，2π］區間內的隨機相位。

微振動模型的優點在于可以根據公式積分推導出加速度模型所對應的理論位移，便于對仿真信號積分處理后的結果進行驗證。

參考微振動信號特性［15］和其他現有微振動模型參數［16］，微振動模型參數體現如下：各階諧波級次呈倍數遞增，一般低階諧波級次振幅較大，是微振動信號中的主導；高階諧波級次振幅較小，通?？梢哉J為前幾階諧波級次生成的信號足以反映微振動信號特性。因此，選用3 階諧波級次用于生成仿真信號，設定變化因子αij和βij在［0.75，1.25］區間內隨機，其他微振動模型信號生成參數如表1所示。

表1 微振動模型信號生成參數Tab.1 Signal generation parameters of micro-vibration model

根據式（15）和表1，生成不同飛輪轉速下的微振動加速度仿真信號。模擬飛輪轉速為2 kr/min 時的微振動加速度仿真信號如圖2 所示，其中：a1，a2，a3為仿真信號各成分；asum為總的仿真信號。

圖2 飛輪轉速為2 kr/min 時的微振動加速度仿真信號（1~1.5 s）Fig.2 Simulation signal of micro-vibration acceleration at flywheel speed of 2 kr/min (1~1.5 s)

2.2 頻譜特征獲取的仿真驗證

以2.1節生成的模擬飛輪轉速為2 kr/min時的微振動加速度仿真信號為對象，加入信噪比為10 dB 的高斯白噪聲，驗證本研究SGHT 方法的頻譜特征獲取能力，同時選取HHT 方法進行對比。為抑制HHT 和SGHT 方法處理后端點效應的影響，HHT和SGHT 方法的輸入均經過極值對稱延拓方法［17］預處理。

使用HHT 方法對預處理信號的處理結果如圖3 所示。

圖3 HHT 方法對預處理信號的處理結果Fig.3 Processing results of simulation signal by HHT

HHT 方法的主要思路是：先對輸入信號進行EMD 分解，再對分解得到的本征模函數（intrinstic mode function， 簡稱IMF）分量做希爾伯特轉換，得到信號分量的瞬時頻率。因為本研究更加關注頻域上微振動能量的整體分布特征，所以選用希爾伯特邊際譜。

由圖3 可以看出，HHT 方法難以對仿真信號進行分解，其邊際譜中可辨別的頻率峰值出現在了11，19，23，58 和63 Hz，其中只有23 Hz 的信號分量與仿真信號某成分特征相符，其他信號分量均不同程度出現了模態混淆現象，得到的頻譜特征結果無法反映輸入信號的基本特征。

使用SGHT 方法對預處理信號的處理結果如圖4 所示。

圖4 SGHT 方法對預處理信號的處理結果Fig.4 Processing results of simulation signal by SGHT

由圖4 可以看出，SGHT 方法能夠準確地對輸入信號進行分離、降噪及識別，得到特征頻率分別為26，58 和84 Hz 的信號分量，同時分離得到的各SGC 分量與混合前的各信號成分基本相符，可反映原始信號成分的物理意義。與HHT 方法相比，SGHT 方法在降噪和分解方面具有顯著優勢，能夠有效處理微振動加速度仿真信號，因此SGHT 算法在頻域特征獲取方面符合微振動信號處理要求。

2.3 時域信號處理的仿真驗證

同樣以2.1節生成的模擬飛輪轉速為2 kr/min 時的微振動加速度仿真信號為對象，加入信噪比為10 dB 的高斯白噪聲，驗證本研究SGHT 方法的時域信號處理能力，即根據加速度仿真信號處理得到線位移估計，并與工程常用處理方法進行對比。

用于對比的常見處理方法包括：加速度處理方法高通濾波；積分方法辛普森公式積分； 位移處理方法多項式去趨勢和EMD。組合這3 種方法得到整體信號處理方法，同時還有整體處理信號的小波積分方法。高通濾波設置低頻截止頻率為10 Hz。小波積分參考文獻［3］中的參數，設置最低截止頻率為5 Hz，最高截止頻率為500 Hz，窗函數為高斯窗（σ=10）。同時，為抑制EMD 和SGHT 方法處理后端點效應的影響，EMD 和SGHT 的輸入經過極值對稱延拓方法預處理。

對微振動加速度仿真信號采用選定的加速度處理方法、積分方法及位移處理方法進行處理，比較處理得到的線位移估計與微振動模型直接得到的理論位移，各方法得到的線位移估計誤差結果如表2所示。

表2 各方法得到的線位移估計誤差結果Tab.2 The error results of the estimated linear displacement processing by each method

表2 中，平均差值誤差定義為線位移估計xest(t)與理論位移xideal(t)的差值在時程內的峰值相對于理論位移xideal(t)峰值的誤差平均值，即

累計相對誤差定義為線位移估計與理論位移差值絕對值的總和與理論位移絕對值總和的比值，即

其中：xest(i)，xideal(i)分別為xest(t)，xideal(t)時程的位移采樣值；Npoint為信號的時間采樣點數。

部分方法線位移估計與理想位移的絕對誤差如圖5 所示。

圖5 部分方法線位移估計與理想位移的絕對誤差（0~0.1 s）Fig.5 The absolute error between the estimated linear displacement and ideal displacement of partial methods（0~0.1 s）

由表2 和圖5 可以看出：未處理的加速度信號在二次積分后，即便以多項式去趨勢方法處理位移，得到的結果依然出現了較大的誤差偏差，這是因為二次積分會放大未處理信號低頻段的微小噪聲，導致積分結果無法體現實際位移；高通濾波和小波積分針對這一問題，通過對低頻段加速度信號的過濾和衰減，使積分誤差降低至一定范圍，但缺點是會損失低頻段信號的部分有效信息；EMD 方法通過分解分離信號中的長周期趨勢項，相較于多項式去趨勢方法的簡單擬合具有一定優勢，但EMD 方法的處理結果受端點效應影響，造成兩端的絕對誤差較大；將SGHT 方法同時應用于加速度處理和位移處理，其能夠在不損失低頻段信號有效信息的基礎上，對加速度信號中的噪聲進行有效的處理，得到與理論誤差最為接近的線位移估計結果。因此， SGHT 方法在時域信號處理方面符合微振動信號處理要求，可以有效得到線位移估計數據，且優于常用的其他方法。

3 工程應用

使用SGHT 方法對ASO-S 衛星/ FMG 載荷地面微振動信號進行處理，結果如圖6 所示。

圖6 SGHT 方法對ASO-S/FMG 微振動信號的處理結果Fig.6 Processing results of ASO-S/FMG micro-vibration signal by SGHT

根據ASO-S 衛星微振動試驗大綱，在FMG 載荷安裝面上共布置4 個三軸線加速度傳感器，負責記錄飛輪不同轉速工況下的微振動試驗數據。在微振動試驗結束后，使用SGHT 方法對各工況下4 個加速度傳感器的3 個軸向加速度數據進行頻域處理。其中四飛輪在2 kr/min 穩定轉速工況下1 個加速度傳感器的z方向經處理后的數據見圖6。由圖可以看出，SGHT 方法能夠對實際微振動加速度信號進行有效的分離、降噪及識別，識別到特征頻率分別為59，65 和345 Hz 的信號分量，其中能量較大的65 Hz信號分量的頻率約為該工況飛輪轉速的2 倍，是該測點的主要微振動來源。根據這些頻譜信息，可以進一步指導載荷的結構設計，以避免在相應頻率附近發生共振。

同時，為了判斷平臺擾動對光學成像的具體影響，需要根據試驗測得的微振動加速度信號估計載荷安裝面的角位移，即進行時域處理。首先，用SGHT 作為加速度處理方法，對各加速度傳感器各軸向加速度數據進行處理；其次，用辛普森公式作為積分方法，對處理后的加速度信號進行二次積分；然后，再次用SGHT 方法作為位移處理方法，分離趨勢項，得到各測點的微振動線位移估計；最后，用最小二乘法擬合各測點形成平面，得到測點安裝面即載荷安裝面的角位移估計。對四飛輪同時工作在2 kr/min 穩定轉速工況的數據進行處理，經SGHT 方法處理得到的載荷安裝面微振動角位移估計如圖7 所示。由圖可以看出，經SGHT 方法處理得到的該工況下載荷安裝面微振動角位移估計均在0.015"以內，進一步分析角位移是以10 Hz 特征頻率為主的周期信號，這一角位移估計可以作為載荷外部的平臺擾動數據，進一步驗證FMG 載荷穩像系統的擺鏡控制算法，判斷微振動平臺擾動對載荷成像是否構成影響。

圖7 經SGHT 方法處理得到的載荷安裝面微振動角位移估計Fig.7 The estimated angular displacement of micro-vibration of payload surface processed by SGHT

通過對不同轉速工況下的試驗數據進行處理，所得到的角位移估計還可以用于確定載荷在軌工作時合適的飛輪轉速工況，減小工況對載荷性能的影響。

4 結 論

1） 在頻域特征獲取方面，SGHT 方法可以在有噪聲干擾的情況下，分解、重構信號并輸出各信號分量及對應的頻譜特征。

2） 在時域信息處理方面，使用SGHT 方法作為加速度和位移處理方法，能夠得到較為準確的線位移估計結果。

3） 使 用SGHT 方法 處 理ASO-S 衛 星FMG 載荷的地面微振動試驗信號，處理得到的頻譜特征和角位移估計可以進一步指導載荷的研制工作。

4） SGHT 方法作為一種新興信號處理方法，目前還存在計算規模大、易受端點效應影響等缺點，有待進一步研究和完善。