龔明程 張學(xué)薇 張樹楨

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

0 引言

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

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

1 基于信號(hào)整數(shù)抽取的頻譜細(xì)化算法簡(jiǎn)介與原理

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

對(duì)z（n）進(jì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為長(zhǎng)度沿坐標(biāo)軸ωx向右進(jìn)行平移，然后再將每次平移后頻譜函數(shù)的幅值相加起來乘以1/D而得出的。再?gòu)摩貀的角度觀察，Y（ejωy）可以理解為Y（ejDωx）順頻率軸擴(kuò)展D倍。以D=3 為例，上述整個(gè)過程如圖1 所示，圖1（a）～圖1（c）為產(chǎn)生混疊的情況，圖1（d）～圖1（f）為不產(chǎn)生混疊的情況，因此，抽取前必須加上防混疊低通濾波器。

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

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

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

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

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

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

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

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

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

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

圖2 Butterworth 低通濾波器幅值響應(yīng)（D=10）

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

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

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

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

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

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

1.3 算法的程序步驟

1.3.1 輸入控制參數(shù)

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

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

1.3.2 確定細(xì)化參數(shù)

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

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

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

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

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

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

1.3.3 時(shí)域復(fù)調(diào)制

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

1.3.4 降重抽取

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

1.3.5 頻域調(diào)整

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

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

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

細(xì)化再迭代。如果細(xì)化后的頻譜分辨率仍不理想，重就復(fù)步驟2、3、4 和5，直至達(dá)到目標(biāo)頻譜精度。

2 仿真算例

以下2 個(gè)算例均在8 核16 線程CPU、3.2 GHz 主頻和16 G內(nèi)存的計(jì)算機(jī)上進(jìn)行計(jì)算。

2.1 含多個(gè)密集頻率成分的頻譜細(xì)化分析

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

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

表1 仿真信號(hào)6 次細(xì)化參數(shù)和頻譜

圖3 仿真信號(hào)的細(xì)化頻譜

2.2 直升機(jī)振動(dòng)信號(hào)的頻譜細(xì)化分析

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

圖4 駕駛員座椅地板振動(dòng)信號(hào)和旋翼轉(zhuǎn)速信號(hào)

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

圖5 駕駛員座椅細(xì)化頻譜和旋翼轉(zhuǎn)速概率密度分布

3 結(jié)語

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

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