基于頻率成比例重分配線性Chirplet變換的時頻分析方法*

2023-12-20 12:38:40岳子毫王征兵李志遠雷歡歡

機電工程 2023年12期

岳子毫,裴幫,王征兵,李志遠,雷歡歡

(鄭州機械研究所有限公司,河南鄭州 450001)

0 引言

齒輪箱是傳動系統中的一個重要部件,廣泛應用于各個領域,其性能直接影響整個設備的健康狀態。因此,齒輪箱的狀態監測和故障研究對保證機械設備的安全穩定運行具有重要意義[1]。

目前,振動分析是齒輪箱狀態檢測與故障診斷中應用最廣泛和最有效的方法[2]。齒輪箱在運行時受到轉速和負載變換的影響,因此,其所收集到的振動信號是非平穩的,且信號分量的振幅和頻率隨時間變換。時頻分析方法是目前處理非平穩信號的最有效的方法之一,其已經在工程上得到了廣泛的應用[3-5]。

時頻分析可以將一維時域信號映射至二維時頻域中,再從時頻域中提取非平穩信號的時變幅度和頻率特性。傳統的時頻分析方法(如短時傅里葉變換和小波變換)是利用計算信號與時頻基之間的內積,從而得到信號中的瞬時特征[6];但是受到海森堡不確定原理的影響,其所得到的時域分辨率和頻域分辨率不能同時達到最佳。Wigner-Ville分布采用雙線性積分變換,其可以獲得更高的時頻域聯合分辨率;但是存在交叉項的干擾問題[7]。

傳統的各類時頻分析方法存在各自的限制問題。因此,為獲得更好的時頻域表示,近些年許多學者在此方面進行了研究,并提出了大量有效的方法。

根據內積定理,基函數與信號的非平穩特征越匹配,得到的時頻表示集中度也越高。MANN S等人[8]提出了Chirplet變換,使用參數Chirprate獲得線性Chirplet基,以此來匹配信號的瞬時頻率特征,相較于傳統的時頻分析方法,其能夠獲得更好的時頻表示;但是Chirplet變換的線性調頻基只能匹配線性瞬時頻率。由Weierstrass逼近定理可知,非線性函數可以被足夠階的多項式近似。為了處理非線性信號,PENG Zhi-ke的團隊[9-10]提出了多項式Chirplet變換,通過構造多項式基來匹配信號的非線性瞬時頻率;但是其在處理多分量信號時,不能使多項式基與所有分量匹配,且多項式階數直接影響計算的復雜度,需要權衡匹配精度和計算復雜度的關系。為了解決這個問題,YANG Yang等人[11-12]又提出了樣條Chirplet變換、廣義鶯變換;但是這些方法不能同時匹配所有分量的瞬時頻率。

YU Gang等人[13]提出了一般線性Chirplet變換(general linear Chirplet transform, GLCT)方法,通過融合多個Chirplet變換產生的結果中能量聚集度的部分來獲得時頻表示,所得到結果的能量聚集度高于任何單個Chirplet變換的結果;但是GLCT在處理瞬時頻率軌跡間隔較小的信號時,分離能力較差[14]。

以上方法都采用了固定的窗口長度,從而導致只能得到固定的分辨率。GUAN Peng-fei等人[15]提出速度同步線性Chirplet變換(velocity synchronous linear Chirplet transform,VSLCT),即通過構造一組與軸速度同步的線性Chirplet基函數,以此來匹配瞬時頻率的軌跡,構造了非正交基,并考慮了時變窗口,其可以清楚地表征所有軸速同步頻率分量;但是該方法計算復雜度大,且時頻表示會出現中斷現象。

時頻方法在處理信號所得到的時頻域表示時,存在拖尾現象。因此,為了提高時頻聚集度,一些學者提出了一些后處理方法。

AUGER F等人[16]提出了重新分配方法(reassignm-ent method, RM),將時頻分布從原始位置重新分配到局部范圍譜的質心,從而提高時頻聚集度;然而RM的過程是不可逆的,不能實現信號重建目的。在RM的基礎上,DAUBECHIES I等人[17]提出了同步擠壓小波變換(synchronous squeezing wavelet transform, SST),其可以較好地完成信號的重建工作。隨后,OBELIN T等人[18]提出了SST的傅里葉變換版本(Fourier-based synchro squeezing transform,FSST),在重新分配過程中,SST會將噪聲引入時頻結果。YU Gang等人[19]提出了同步提取變換(synchronous extraction transfor-mation,SET),利用狄拉克函數和估算頻率構建了同步提取算子,只提取與瞬時頻率最相關的成分,避免了噪聲的干擾,相較于SST,其能獲得更高的時頻聚集度。

然而,以上這些方法在處理強時變信號和高噪聲信號時的效果并不理想。

為了對齒輪進行故障診斷,同時針對上述問題,筆者從齒輪故障振動信號的調制現象出發,討論基函數對時變信號分析的影響,提出一種時頻方法-頻率成比例重分配Chirplet變換(PFSRLCT)方法,并分別在模擬信號和實際齒輪故障信號上對PFSRLCT方法進行驗證。

1 齒輪故障信號調制現象

齒輪在正常工況下嚙合時,其振動信號由齒輪的嚙合頻率及其諧波組成,即:

(1)

式中:Ai為i階諧波的幅值;fm為齒輪嚙合頻率;φi為i階諧波的相位。

齒輪嚙合傳動中,負載、剛度和轉速的變化以及故障的發生,都會使得齒輪振動信號產生調制現象。由此,齒輪的振動信號可以表示為:

(2)

式中:ai(t)為調幅函數;bi(t)為調頻函數。

齒輪故障信號產生調幅調頻現象時,時域和頻域如圖1所示。

圖1 齒輪故障信號的調幅調頻現象Fig.1 Amplitude modulation and frequency modulation of gear fault signal

對于齒輪振動信號,嚙合頻率為載波頻率fm,齒輪的旋轉頻率為調制頻率fr。當齒輪發生斷齒等局部故障時,振動信號會被調幅,在頻譜上產生以載波頻率與調制頻率的和頻及差頻。如果調制信號不是單一信號,而是具有諧波成分的信號,則會在載波頻率及其諧波兩側出現以調制頻率為間隔的一系列邊頻帶。局部故障還會產生調頻現象,在頻譜上會產生以載波頻率為中心,以調制頻率為間隔,對稱分布的無限多對的邊帶。在實際的齒輪系統中,調幅和調頻一般同時存在。

調幅調頻現象同時存在時,振動信號的頻譜如圖2所示。

圖2 調幅調頻同時存在的信號頻譜Fig.2 Signal frequency spectrum when AM and FM exist at the same time

由于邊頻成分相位的不同,矢量相加后會在頻譜上形成復雜的不對稱調制邊帶。同時,齒輪的旋轉頻率與載波頻率存在以下關系:

fr=fm/Z

(3)

式中:Z為齒輪齒數。

由上述分析可知,齒輪振動信號中各個頻率分量的頻率成比例。

2 Chirplet變換

信號x(u)òL2(R)的Chirplet變換可以表示為:

(4)

式中:A(u)為信號x(u)的Hilbert變換,表示為s(u)=x(u)+jH[x(u)];g(t)為窗口函數,通常采用高斯函數g(t)=exp(-t2/(2σ2)),其中σ為高斯窗口參數;φ(f,u,tc)為相位函數。

φ(f,u,tc)可表示為:

φ(f,u,tc)=C(u-tc)2/2

(5)

式中:C為Chirprate。

Chirplet變換本質上是一種加窗變換,在時間窗口內通過旋轉Chirplet,使得其與瞬時頻率軌跡匹配,以此來獲得更高的時頻能量聚集度。

與短時傅里葉變換(short time Fourier transform,STFT)相比,其不同在于Chirplet變換利用Chirplet基與瞬時頻率軌跡匹配,而不是時頻基。時頻能量聚集度和Chirplet與瞬時頻率軌跡重合程度成正比,當Chirplet與瞬時頻率軌跡完全重合時,應該存在Chirplet旋轉角θ的正切值等于目標瞬時頻率的軌跡的斜率。

Chirplet的旋轉角度可以表示為:

φ′(f,u,tc)=dφ/du=C(u-tc)

(6)

-tan(θ)=dφ′/du=C

(7)

觀察式(4)可以得出,當C=0,Chirplet變換退化為短時傅里葉變換。

Chirplet變換相較于短時傅里葉變換能量聚集度更高的原因如圖3所示。

圖3 Chirprate對時頻聚集度的影響Fig.3 The influence of CR on TF energy concentration degree

由圖3可知:短時傅里葉變換的時頻基與瞬時頻率軌跡僅在時間窗口寬度中心處匹配,在窗口的其他位置與真實時頻特征不匹配,能量分布較為分散。

筆者對單分量非平穩信號進行Chirplet變化時,Chirplet基瞬時頻率軌跡匹配如圖4所示。

圖4 單分量信號的Chirplet變換Fig.4 Chirplet transform of single component signal

由圖4可知:當所考慮的信號包括時變頻率分量時,任何單個Chirplet變換都不能與瞬時頻率軌跡完美匹配。

針對上述問題,YU Gang提出了GLCT,即在每個時間窗口內尋找與該時間窗口內瞬時頻率軌跡匹配的Chirprate,結合多個Chirplet變換的結果,引入參數α,其公式如下:

(8)

式(4)改寫為:

(9)

當N=1時,GLCT退化為Chirplet變換。

GLCT有效地解決了瞬時頻率軌跡斜率隨時間變換的問題。筆者采用時變Chirplet匹配頻率軌跡,在處理單分量信號時可以獲得良好的能量聚集度。

GLCT處理多分量非平穩信號如圖5所示。

圖5 多分量信號的GLCTFig.5 GLCT transform of multi-component signal

由圖5可知:當處理多分量信號時,如齒輪箱瞬時頻率成比例的多分量信號,在每個時間窗內的各個分量信號的瞬時頻率具有不同的斜率,GLCT不能完全匹配所有的分量信號。此外,GLCT不適合處理頻率間隔較近的多分量信號。

3 頻率成比例重分配線性Chirplet變換

為了解決時頻分析存在的問題,筆者提出了頻率成比例重分配線性Chirplet變換(PFSRLCT)方法。

首先,替換Chirplet變換的核函數,得到成比例頻率線性Chirplet變換(proportional frequency linear Chirplet transform,PFLCT);然后,提取時頻脊,獲得同步重分配算子(synchronous reassignment operator,SRO),對時頻結果進行了同步分配,進一步提高了能量集中度。

3.1 成比例頻率線性Chirplet變換

修改式的相位函數為:

φ(f,u,tc)=fC(u-tc)2/2

(10)

式中:f為頻率。

由此,式(4)改寫成:

(11)

由Ville理論[20]和泰勒公式可得:

(12)

f(u)=f(tc)+f′(tc)(u-tc)

(13)

將式(12)和式(13)代入式(11),可得:

|PFLCT(f,tc)|=

g(u)exp(-j2πf(tc)Cu2)exp(-jωu)du|=

exp(-j(ω-2πf(tc))u)du|

(14)

(15)

由上式可知,當C滿足C=f′(tc)/2f(tc)時,ω在2πf(tc)處|PFLCT(f,tc)|取得最大值,即獲得最高能量濃度。由此可以證明,PFLCT可以處理具有瞬時頻率成比例的信號。

另外,PFLCT還可用于處理具有瞬時頻率成比例的多分量信號,證明過程如下:

假設多分量信號f(t)=f1(t)+f2(t)且?f1(t)=af2(t),那么對于信號f(t)進行PFLCT變化可得:

(16)

(17)

由式(17)可證PFLCT可用于分析具有成比例瞬時頻率的多分量信號。

(18)

式中:N為從-π/2到π/2的整個角度跨度中的段數。

即對于?tc,?tanαk=f′(tc)/2f(tc),k∈[1,N],當時間窗口內時頻基與瞬時頻率軌跡最匹配時,時頻域的能量集中在各分量的對應頻率上,這時時頻域的峭值達到最大值。

因此,可以使用峭值作為指標來確定C的值:

(19)

式(18)中,N的值越大,式(19)中C的值越精確,時頻基和瞬時頻率軌跡匹配程度越良好;但是過大的N值會增加計算負擔。

筆者在該項研究中設置N為50。

PFLCT的算法流程如下:

算法1:PFLCT輸入A(u),σ,fs,Nfor i=1:NTFRs(t,f,i)=∫+∞-∞ A(u)g(u-t)exp(-jω(u-t))exp(-j2πftan-π2+iN+1π()(u-t)2)duend forCI=Kurtosis(TFRs,t)Index=argmaxc(CI)TFR=TFRs(index)輸出TFR,index,t,f

3.2 同步重分配

對于多分量信號,其中每個分量信號的瞬時頻率應該有足夠的間隔,才能從時頻域中提取出各分量時頻脊線,即:

|fi(t)-fi+1(t)|>2Δ

(20)

式中:fi(t)為第i個分量信號的頻率,i∈Z;Δ為固定頻率間隔。

目前,已經有了大量的研究來解決時頻脊線提取的問題。THAKUR G等人[21]開發了與式(20)的分辨率相關的另一種簡單而有效的方法,其表達式如下:

(21)

式中:(t,φk(t))為時頻域內的估計瞬時頻率軌跡;λ,β為曲線平滑度和能量最大化的權衡參數。

該方法采用貪婪算法尋找局部最優解,其具體步驟為:

1)將時頻信號按時間方向分割為若干片段,并確定每個片段的搜索起點;

2)以最大化式為衡量標準,使用向前向后的方式搜索瞬時頻率脊線;

3)每獲取完整的一條瞬時頻率脊線后,在原信號去除該分量;

4)重復步驟1)～步驟3),直到提取出所有瞬時頻率脊線。

尋找時頻脊線的算法流程如下:

算法2:尋找時頻脊線輸入TFR,N,β,M,N[flen tlen]=size(TFR)for i=1:Mfor k=1:NIF(tk)=argmaxf(TFR(k·tlen/(N+1),f))IF(tk-1)=argmaxf(TFR(k·tlen/(N+1)-1,f))for j=k+1:tlenIF(tj)=argmaxf∈[IF(tk-1)-Δ,IF(tk-1)+Δ](|TFR(f,tj)|2-β(f-2IF(tj-1)+IF(tj-2))2)end forfor j=0:k+1IF(tj)=argmaxf∈[IF(tk+1)-Δ,IF(tk+1)+Δ](|TFR(f,tj)|2-β(f-2IF(tj+1)+IF(tj+2))2)end forend forIFs(:,M)=IFTFR(t,(IF(t)-Δ,IF(t)+Δ))=0clear IFend for輸出IFs

筆者利用提取到的各分量的時頻脊線和狄拉克函數,使得時頻能量僅在瞬時頻率軌跡處分布。

因此,可得到同步分配算子SRO如下:

(22)

式中:IF(t)為時頻脊線。

PFSRLCT可以表示為:

PFSRLCT(f,t)=PFLCT(f,t)·SRO(f,t)

(23)

重分配過程是僅將每條時頻脊線上的局部最大時頻幅值重新分配給新的時頻域表示。因此,重分配后的時頻域表示具有更高的時頻聚集度。

筆者對時頻結果進行重分配的算法流程如下:

算法3:重分配TFR輸入TFR,IF,M,Δkk=1:hlength:hlength·llengthindex=kk+IF-1;SRO=zeros(size(TFR))for i=1:Mindex1=repmat(index,Δ+1,1)+linspace-Δ2,Δ2,Δ+1()TSRO(index1)=1end forTFR=TFR.?SRO輸出TFR

4 模擬信號驗證

筆者使用模擬信號來驗證PFSRLCT的有效性,并將其與其他方法進行了對比。

非線性多分量模擬信號模型如下:

(24)

該多分量信號包括3個瞬時頻率成比例的非線性信號,采樣頻率為1 000 Hz,采樣時長為3.5 s,筆者設置信噪比SNR=-1 dB。

該多分量信號的信號波形如圖6所示。

圖6 模擬信號的時域波形Fig.6 Time domain waveform of analog signal

4.1 表現能力驗證

筆者利用式(24)建立了模擬信號的分析結果,如圖7所示。

圖7 模擬信號分析結果Fig.7 Analysis results of analog signals

圖7(a)為信號的瞬時頻率軌跡,圖7(b)為應用PFLCT獲得的時頻域結果。對比兩者結果可以發現:PFLCT得到的時頻表示可以清楚地分辨出各個分量信號,證明該方法對于處理間隔緊密的瞬時頻率成比例多分量信號具有良好的效果。

圖7(c)為對PFLCT所得時頻結果進行脊提取得到的各分量的瞬時頻率,可以看出:圖中所提取的瞬時頻率軌跡與真實瞬時頻率基本吻合;筆者利用所提取的時頻脊線獲得同步重分配算子SRO,對PFLCT進行重分配后所得時頻域表示如圖7(d)所示。由圖7(b)中的時頻結果比較可知,重分配后的時頻結果能量更加集中。

為了評估PFSRLCT方法的表現能力,筆者采用了一些傳統方法和先進方法STFT、GLCT、SET、VLSCT、SST、SST2進行了比較。其中,STFT、GLCT、SET、SST、SST2的窗寬均設置為0.5 s,VLSCT的Nwl設置為30。

所得結果如圖8所示。

通過觀察圖8可知:STFT受到拖尾效應的影響,所得時頻域表示較為模糊,時頻集中度較低;GLCT在處理f1、f2間隔緊密的頻率時,發生了頻率交叉,難以分辨出頻率軌跡;SET結果時頻集中度高,但f3軌跡在2 s后發生了模糊;VLSCT的時頻表現良好,但是在2.5 s和3 s之間發生了中斷。

圖8 各種時頻方法的比較Fig.8 Comparison of various TFA methods

圖8(e)、圖8(f)分別為SST和SST2的結果,這兩種方法相對于STFT有較好的時頻集中度,但是在幅值上有所削弱,在頻率間隔較近時仍有模糊。

為了驗證PFSRLCT的抗噪性能,筆者引入了Rényi熵來衡量時頻域的集中度。

Rényi熵由下式確定:

(25)

其中:α設置為3。

Rényi熵可以量化TFR的清晰度,Rényi熵越高意味著時頻表示越模糊。

通過式(25)可計算得到各種時頻方法的Rényi熵,如表1所示。

表1 各種方法的Rényi熵對比

由表1可知:GLCT得到時頻結果的Rényi熵最大,為10.393 7;PFSRLCT所獲得的Rényi熵最小,為4.745 1,并且明顯低于其他時頻方法。

由此可以證明,PFSRLCT所獲得時頻能量集中度最高,相對于其他時頻方法有明顯優勢。

4.2 噪聲魯棒性

為了檢測該方法對于噪聲的魯棒性,筆者將具有不同信噪比的高斯噪聲添加到信號中,使用頻率估計誤差來量化頻率估計的準確性。

頻率估計誤差表示如下:

(26)

在提取時頻脊的過程中,筆者提取7條時頻脊線,并通過人工選取與真實瞬時頻率最符合的3條脊線;利用式(26),得到模擬信號在信噪比為-2 dB～10 dB時頻率估計誤差,并將其與STFT、SST、VSLCT結果進行對比。

4種方法的頻率估計誤差,如圖9所示。

圖9 頻率估計誤差對比Fig.9 Comparison of frequency estimation errors

從圖9中可以看出:在信噪比較高時,SST2表現效果不佳,在信噪比高于1 dB時,表現效果與VSLCT相似,次于PFSRLCT的誤差;由于STFT時頻結果的能量局部帶寬更大,干擾了脊提取算法,其結果的頻率誤差相對于其他3種算法最大,PFSRLCT所得到的頻率估計誤差最小,同時受信噪比的影響小。

在信噪比為-1 dB時,筆者采用4種方法得到了真實頻率和估計頻率,如圖10所示。

圖10 真實瞬時頻率軌跡與估計瞬時頻率軌跡的比較Fig.10 Comparison between real IF trajectory and estimated IF trajectory

由圖10可以看出:PFSRLCT所得到的估計頻率與真實頻率符合度最高;VSLCT所得到的估計頻率在2 s左右頻率間隔較小處和2.5 s處有較大突變;STFT與SST2所得到的估計頻率在2 s到2.5 s之間頻率間隔較近處有較大誤差,STFT所估計的f3受能量局部帶寬影響較大,所得到的估計頻率在真實頻率上下浮動;受到噪聲的影響,SST2所估計的f3在2.5 s后有較大的偏離。