基于信號整數(shù)抽取的頻譜細化算法在直升機振動信號分析中的應用

龔明程 張學薇 張樹楨

（航空工業(yè)直升機設計研究所，天津 300000）

0 引言

在旋轉機械健康監(jiān)測和故障診斷中，需要關注的頻譜分量往往集中在一個狹窄的頻帶內，例如拍振信號就是一種典型的密集頻率信號[1-2]。通常情況下，直升機的旋翼轉速是比較穩(wěn)定的，也可以看成集中在窄帶的頻率。由于受計算機內存和計算能力的限制，因此現(xiàn)實中會選擇1 段局部信號進行FFT 變換，該方法以頻譜分辨率為代價，減少了頻譜分析中所需信號的樣本數(shù)量，提高了分析效率。因此，既能獲得高質量頻譜分辨率，又比完整信號FFT 變換耗時少、效率高且計算量小的信號分析方法在工程中是必要的。

針對密集頻率成分的識別問題，專家們已經提出多種頻譜細化和校正方法，常見的方法有線性調頻Z 變換法[3]、FFT+FS 法[4]和復調制頻譜細化法（ZOOM-FFT）等。其中，ZOOM-FFT 法具有計算精度高和速度快的優(yōu)點，適用于FFT變換點數(shù)不多的情況[5-6]?；谝陨纤悸罚撐奶岢隽艘环N高效的區(qū)分信號密集頻率成分的算法，并從數(shù)學上完整證明了該方法的原理。在同等處理器計算能力瓶頸的限制下，該方法的頻譜質量比普通FFT 更高，得到工程中所關心頻率附近的高精度頻譜。

1 基于信號整數(shù)抽取的頻譜細化算法簡介與原理

1.1 信號整數(shù)抽取的頻域特性

設有1 個無窮離散序列x（n），其中n∈Z，則用正整數(shù)D抽取該序列在時域上的表現(xiàn)為y（n）中每隔D-1 個抽樣點抽取1 個樣點，然后按次序組成序列x（n），該過程可用數(shù)學表達式描述，如公式（1）所示。

構造1 個序列z（n），該序列由脈沖序列和序列x（n）相乘而來，如公式（2）所示。

式中：i為整數(shù)；n為整數(shù)。

顯然，序列x（n）、y（n）和z（n）在時域上的關系如公式（3）所示。

設x（n）、y（n）和z（n）的抽樣頻率分別為fx、fy和fz，可以得到公式（4）。

對各自抽樣頻率歸一化的數(shù)字頻率變量ωx、ωy和ωz來說，可定義如公式（5）所示。

則3 個數(shù)字頻率變量由公式（6）關聯(lián)。

對y（n）做離散時間傅里葉變換（DTFT），如公式（7）所示。

式中：j為虛數(shù)單位；z（m）為構造的序列；Z為進行離散時間傅里葉變換后的頻譜函數(shù)。

由文獻[7]和文獻[8]可知，沖激函數(shù)具有尺度變換和取樣特性，如公式（8）所示。

式中：a為非零實數(shù)；δ為沖擊函數(shù)；q（t）為一個在t=t0處連續(xù)的普通函數(shù)。

利用公式（8），結合只有當i=0 時，ω-2πi才能落在[-π，π]的積分區(qū)間的特性，可得到恒等式，如公式（9）所示。

c（n）又可以為公式（11）。

根據(jù)離散時間傅里葉變換成對的性質，又可得到恒等式，如公式（12）所示。

另一方面，直接對c（n）進行DFTF 變換，如公式（13）所示。

聯(lián)立公式（12）和公式（13），可以得到公式（14）。

對z（n）進行DTFT 變換，如公式（16）所示。

式中：k為正整數(shù)。

將公式（16）和公式（6）代入公式（7），得到公式（17）。

式中：y（ejωy）為頻譜函數(shù)。

先從ωx的角度觀察公式（17），Y（ejDωx）是先將y（ejωx）依次以0、2π/D、…、2（D-1）π/D為長度沿坐標軸ωx向右進行平移，然后再將每次平移后頻譜函數(shù)的幅值相加起來乘以1/D而得出的。再從ωy的角度觀察，Y（ejωy）可以理解為Y（ejDωx）順頻率軸擴展D倍。以D=3 為例，上述整個過程如圖1 所示，圖1（a）～圖1（c）為產生混疊的情況，圖1（d）～圖1（f）為不產生混疊的情況，因此，抽取前必須加上防混疊低通濾波器。

圖1 初采樣序列和再抽取序列的頻譜（D=3）

由于Y（ejωy）的周期性，因此可以僅考慮單周期ωy∈[-π，π]，則ωx∈[-π/D，π/D]。此時只有k=0，X（ej（ωx-2πk/D））才會有意義，因此也只需要對其進行取k=0 的單周期分析，公式（17）變?yōu)楣剑?8）。

離散時間序列DTFT 變換后的結果是關于頻率的連續(xù)函數(shù)，無法使用計算機處理，必須利用DFT 變換對其頻域進行離散化，由于工程實際中的信號均為有限長信號，因此可以設采樣后的有限序列為x1（lx），lx=0，1，2，…，N，即x1（lx）為無窮離散序列x（n）的一部分，則x1（lx）用正整數(shù)D抽取后的序列y1（ly）為公式（19）。

根據(jù)文獻[9]可知，序列x1（lx）的 DFT 和x（n）的DTFT變換之間的關系、序列y1（ly）的 DFT 和y（n）的DTFT 變換之間的關系如公式（20）所示。

式中：RND（kx）和RN（ky）為矩形序列。

將公式（20）代入公式（18），并進行如公式（21）所示的變形。

由公式（21）可知，Y1（ky）的N條譜線和X1（kx）中的前N條譜線是完全相等的。因此，可以通過抽取后的信號y（ly）去間接求解原始信號的頻譜，并且抽取后的信號的長度縮短，以較少的DFT 變換點數(shù)得到了原始信號的頻率分辨率。

值得注意的是，為了保證降重抽取后的信號不產生混疊失真，需要在抽取前加上防混疊低通濾波器，由文獻[10]和文獻[11]可知，防混疊低通濾波器的理想頻率響應Hd（eiω）須滿足公式（22）。

該文提供了2 種抽取數(shù)字濾波器設計方法，第一種方法基于Butterworth 生成IIR 低通數(shù)字濾波器，通帶最大衰減As、阻帶最小衰減Rp、歸一化通帶頻率和歸一化阻帶頻率4 個參數(shù)的取值如公式（23）所示。

根據(jù)Butterworth 濾波器的幅度平方函數(shù)可以求解濾波器階次L和歸一化截止頻率，計算方法如公式（24）所示。當D=10 時，第一種方法對應濾波器的幅值響應曲線如圖2 所示。

圖2 Butterworth 低通濾波器幅值響應（D=10）

第二種方法是直接利用MATLAB 自帶的resample 函數(shù)進行抽取，該函數(shù)內置FIR 抗混疊低通濾波器，并補償濾波器帶來的延遲。因此，使用resample 函數(shù)可以進一步提高程序的運行速率。

1.2 離散傅里葉變換的頻域循環(huán)移位特性

設有限長度的離散時間序列為p0（n）（其中，n=0，1，…，N-1），采樣頻率為fs，則該序列的N點DFT 如公式（25）所示。

同時，譜線間隔記為?f，實際需要細化分析的頻帶區(qū)間記為[fl，fu]，頻帶中心頻率記為fe。由于細化分析時需要將中心頻率fe移頻至頻率軸的零點，不妨先假設頻譜P0（k）向左循環(huán)位移M個樣點后，中心頻率剛好移頻至頻率軸的零點，此時的頻譜可記為P（k），如公式（26）所示。式中：（k+M）N為對k+M進行模N運算。

將公式（26）帶入逆離散傅里葉變換IDFT，如公式（27）所示。

因此，中心頻率fe剛好移頻至頻率軸的零點，等價于用e-i2πMn/N對原序列P0（n）進行復調制。參數(shù)M由fe在原頻譜P0（k）中的位置來確定，取值如公式（28）所示。

1.3 算法的程序步驟

1.3.1 輸入控制參數(shù)

信號和分析參數(shù)：原始信號為x（n）；信號采樣頻率為fs；FFT 變換點數(shù)為nfft。

IIR 低通濾波參數(shù)：通帶最大衰減為As；阻帶最小衰減為Rp。

1.3.2 確定細化參數(shù)

計算需要細化分析的頻帶中心頻率fe。對信號x（n）做點數(shù)為nfft的FFT 變換，得到頻率分辨為?f=fs/nfft的頻譜，確定需要細化分析的頻帶區(qū)間[fl，fu]，計算中心頻率fe，如公式（29）所示。

確定最大細化倍數(shù)D。根據(jù)奈奎斯特抽樣定理可知，細化倍數(shù)（整數(shù)抽取因子）D的選取理論上只需要滿足公式（30）。

考慮當設計的低通濾波器在|ω|∈[π/2D，π/D]的過渡帶時，實際頻率響應Hd（eiω）存在衰減，即理論上計算的細化頻帶區(qū)間[fe-fd/2，fe+fd/2]中只有[fe-fd/4，fe+fd/4]區(qū)間的幅值是真實的?？梢酝ㄟ^強化約束細化倍數(shù)的判斷公式（30）解決該問題，公式（30）修正為公式（31）。

公式（30）可以理解為判定2 倍細化倍數(shù)后的理論細化頻帶是否包括需要細化分析的頻帶[fl，fu]。顯然，這種保守做法可以確保需要細化分析頻帶區(qū)間幅值的真實性。為了滿足正整數(shù)抽取的要求，可以根據(jù)公式（32）求解最大細化倍數(shù)。

式中：floor 為向下取整函數(shù)。

設計低通濾波器：計算歸一化通帶頻率ωp和阻帶頻率ωs，基于Butterworth 模擬濾波器生成IIR 低通數(shù)字濾波器。

1.3.3 時域復調制

將中心頻率fe移頻至頻率軸的零點，先計算抽取/調制需要用到的點數(shù)ndec，由于頻域循環(huán)位移等價于對序列x（n）進行時域復調制，因此調制后x（n）變?yōu)閺蛿?shù)信號x1（n），如公式（33）所示。

1.3.4 降重抽取

按照正整數(shù)因子D進行抽取。經過復調制和低通濾波后，信號頻帶變窄，然后再對其進行x1（n）抽取，等價于用較低的采樣率fdec對x1（n）進行重采樣，如公式（34）所示。

1.3.5 頻域調整

FFT 變換。對抽取后的離散序列x2（）進行點數(shù)為nfft的復數(shù)FFT 計算，頻譜記為X2（k）。

幅值修正和翻轉。fftshift 函數(shù)為MATLAB 內置函數(shù)，作用是將0 頻率分量移到頻譜中心，如公式（35）所示。

頻率刻度恢復。重采樣后的頻率分辨率如公式（30）所示，比第一次的分辨率提高了D倍，如公式（36）所示。

細化再迭代。如果細化后的頻譜分辨率仍不理想，重就復步驟2、3、4 和5，直至達到目標頻譜精度。

2 仿真算例

以下2 個算例均在8 核16 線程CPU、3.2 GHz 主頻和16 G內存的計算機上進行計算。

2.1 含多個密集頻率成分的頻譜細化分析

設1 個加速度仿真信號由5 個不同頻率的正弦信號構成，幅值均為1，初始相位均為0，頻率分別為100.047、100.049、100.050、100.051 和100.054，并添加均值為0、方差為1 的高斯白噪聲，信號的采樣頻率為2 000 Hz，單位為m/s2。

仿真算例中的FFT 變換點數(shù)取為2 048，經過6 次頻譜細化迭代后可精確識別5 個非?？拷恼瓗ьl率，最大誤差為0.000 008 Hz，準確率極高，需要細化倍數(shù)和頻帶區(qū)間的設置見表1，每次細化后的頻譜如圖3 所示。如果傳統(tǒng)方法要達到的同樣的精度，就需要進行點數(shù)為2048×50000 的FFT 變換，對計算機性能的要求較高，而該文提出的方法只需要進行點數(shù)為2 048 的FFT 變換，極大地降低了對硬件內存和處理器性能的要求。

表1 仿真信號6 次細化參數(shù)和頻譜

圖3 仿真信號的細化頻譜

2.2 直升機振動信號的頻譜細化分析

某型民用直升機在某次試飛試驗中，500 s～2 500 s 的右駕駛員座椅地板處加速度信號和旋翼轉速Nr 百分比信號如圖4 所示，加速度信號采樣率為1 024 Hz，以重力加速度g為基本單位，Nr 信號采樣率為25 Hz，單位為百分比。已知Nr=100%時的旋翼一階通過轉速頻率為23.65 Hz。

圖4 駕駛員座椅地板振動信號和旋翼轉速信號

該算例利用該文提出的算法準確提取直升機的旋翼一階通過頻率。其中，細化參數(shù)按照如下方法設置，F(xiàn)FT 變換的點數(shù)為1 000，細化倍數(shù)D=2 000，細化頻帶取[23.03，23.28]，則細化后的頻譜如圖5（a）所示，最高峰對應的頻率為23.151 5 Hz。另外，Nr 信號也是離散點，通過統(tǒng)計的方法計算Nr 的每個值出現(xiàn)的次數(shù)，最后再除以500 s～2 500 s 期間Nr 信號的總長度，得到Nr 的每個值在500 s～2 500 s 飛行中出現(xiàn)的概率，也可以理解為Nr 的概率密度分布，如圖5（b）所示，出現(xiàn)次數(shù)最多Nr 對應的頻率為23.148 3 Hz。圖5（a）和圖5（b）中的峰值對應的頻率相差0.003 2 Hz，再次驗證了該文提出的算法具有非常高的精度。

圖5 駕駛員座椅細化頻譜和旋翼轉速概率密度分布

3 結語

該文提出的基于信號抽取性質的頻譜細化算法，一方面降低了需要分析的信號長度，等價于降低對FFT 變換的點數(shù)的要求，降低了對處理器硬件的要求，變相提高了分析效率；另一方面，如果采用之前FFT 變換的點數(shù)，在不改變計算機硬件性能的情況下，可以提高關心頻帶區(qū)間的頻譜分辨率，極高的精度優(yōu)勢是傳統(tǒng)FFT 分析方法無法比擬的。

借鑒第二個算例，該文提出的算法還可以進一步拓展到直升機排故中，先通過傳統(tǒng)FFT 方法分析故障頻率的大概范圍，然后再用該文提出的算法逐步對關心頻帶進行細化分析，也可以為精準定位和進一步分析異常振動源頭提供依據(jù)。