本發(fā)明涉及轉(zhuǎn)子振動(dòng)信號處理與轉(zhuǎn)子模態(tài)參數(shù)辨識領(lǐng)域，尤其涉及運(yùn)行狀態(tài)下的振動(dòng)模態(tài)分析領(lǐng)域。

背景技術(shù)：

振動(dòng)模態(tài)分析是獲取機(jī)械結(jié)構(gòu)動(dòng)態(tài)特性不可或缺的手段，是振動(dòng)控制、結(jié)構(gòu)狀態(tài)監(jiān)測、減震降噪、機(jī)械結(jié)構(gòu)故障診斷、有限元模型修正及確認(rèn)的基礎(chǔ)。目前，主要方法有：有限元法、基于輸入輸出模態(tài)數(shù)據(jù)的傳統(tǒng)試驗(yàn)?zāi)B(tài)分析法和基于僅有輸出數(shù)據(jù)的運(yùn)行狀態(tài)模態(tài)分析法。有限元法對于求解轉(zhuǎn)子和周圍結(jié)構(gòu)一起組成的旋轉(zhuǎn)機(jī)械問題時(shí)有很突出的優(yōu)點(diǎn)。但實(shí)際工程中，由于結(jié)構(gòu)復(fù)雜邊界條件、結(jié)構(gòu)物理參數(shù)和部件連接狀態(tài)等不確定因素的影響，很難建立準(zhǔn)確的有限元模型。

傳統(tǒng)的試驗(yàn)?zāi)B(tài)分析通常在實(shí)驗(yàn)室內(nèi)完成，試驗(yàn)狀態(tài)易于控制，測量信噪比較高。旋轉(zhuǎn)機(jī)械的試驗(yàn)?zāi)B(tài)分析主要在非工作狀態(tài)下進(jìn)行模態(tài)測試，由于缺少陀螺力矩等因素的影響，與其運(yùn)行狀態(tài)下的測試結(jié)果可能存在較大區(qū)別。傳統(tǒng)的運(yùn)行狀態(tài)模態(tài)分析法，一般要求結(jié)構(gòu)處于具有寬頻特征的激振力作用下，而對于旋轉(zhuǎn)機(jī)械，各旋轉(zhuǎn)部件的運(yùn)轉(zhuǎn)將引起與轉(zhuǎn)速密切相關(guān)的諧波分量，給結(jié)構(gòu)模態(tài)參數(shù)識別造成很大困難。

技術(shù)實(shí)現(xiàn)要素：

為了克服現(xiàn)有技術(shù)存在的問題，本發(fā)明實(shí)施例提供了一種基于階次提取的轉(zhuǎn)子運(yùn)行狀態(tài)模態(tài)分析法，能夠避免由旋轉(zhuǎn)激勵(lì)引起的諧波干擾，提高旋轉(zhuǎn)機(jī)械運(yùn)行狀態(tài)下模態(tài)參數(shù)識別的精度。

為達(dá)到上述目的，本發(fā)明的實(shí)施采用如下技術(shù)方案：

第一方面，本發(fā)明實(shí)施例提供一種基于瞬時(shí)頻率估計(jì)的自適應(yīng)vold-kalman濾波階比跟蹤技術(shù)，所述方法用于轉(zhuǎn)子運(yùn)行狀態(tài)振動(dòng)信號處理，所述方法包括：

針對所述轉(zhuǎn)子振動(dòng)信號，利用瞬時(shí)頻率估計(jì)法計(jì)算出轉(zhuǎn)子轉(zhuǎn)速；

根據(jù)所述轉(zhuǎn)速，結(jié)合自適應(yīng)vold-kalman濾波階比跟蹤技術(shù)計(jì)算得到某一階次信號。

第二方面，本發(fā)明實(shí)施例提供一種將自適應(yīng)濾波階比跟蹤技術(shù)與模態(tài)識別算法相結(jié)合的轉(zhuǎn)子模態(tài)分析法，所述方法用于轉(zhuǎn)子的模態(tài)參數(shù)識別，所述轉(zhuǎn)子系統(tǒng)在運(yùn)行狀態(tài)下不方便施加激勵(lì)，在僅有響應(yīng)的情況下，所述方法包括：

針對所述提取的某一階次信號直接進(jìn)行功率譜分析，對所述功率譜進(jìn)行奇異值分解，得到左奇異向量與右奇異向量；

根據(jù)所述左奇異向量與右奇異向量計(jì)算得到一個(gè)增強(qiáng)功率譜，利用所述增強(qiáng)功率譜計(jì)算出阻尼和固有頻率。

本發(fā)明提供的一種基于瞬時(shí)頻率估計(jì)的自適應(yīng)濾波階比跟蹤技術(shù)，與目前其他階比跟蹤算法相比，本實(shí)施例不需要安裝硬件設(shè)備來測量轉(zhuǎn)速信號，大大減少了工作量。此外，本實(shí)施例還提供一種將自適應(yīng)濾波階比跟蹤技術(shù)與頻域空間域分解法相結(jié)合的轉(zhuǎn)子模態(tài)分析法，在不方便施加激勵(lì)或激勵(lì)未知，僅有響應(yīng)的情況下，具有顯著的優(yōu)勢，且工作量少，計(jì)算量小，運(yùn)算速度快。

附圖說明

圖1為本發(fā)明算法識別轉(zhuǎn)轉(zhuǎn)機(jī)械轉(zhuǎn)子的結(jié)構(gòu)模態(tài)參數(shù)的流程圖；

圖2為轉(zhuǎn)子原信號瀑布圖；

圖3為轉(zhuǎn)子原信號與提取2x信號的時(shí)域?qū)Ρ龋?/p>

圖4為2x信號的瀑布圖；

圖5為2x信號的模態(tài)指示曲線；

圖6為本發(fā)明算法識別轉(zhuǎn)子的前三階彎曲模態(tài)。

具體實(shí)施方式

本發(fā)明實(shí)施例提供了一種基于階次提取的旋轉(zhuǎn)機(jī)械轉(zhuǎn)子運(yùn)行狀態(tài)模態(tài)分析法，能夠避免旋轉(zhuǎn)激勵(lì)引起的諧波干擾，從而提高旋轉(zhuǎn)機(jī)械運(yùn)行狀態(tài)下模態(tài)參數(shù)識別的精度。

為達(dá)到上述目的，如圖1所示，本發(fā)明的實(shí)施采用如下步驟：

步驟一：對時(shí)域信號進(jìn)行時(shí)頻分析，畫出瀑布圖；

步驟二：由瞬時(shí)頻率估計(jì)法估計(jì)轉(zhuǎn)軸轉(zhuǎn)速；

步驟三：自適應(yīng)vold-kalman濾波階比跟蹤算法對某一階次信號進(jìn)行階次提取；

步驟四：對提取的階次信號使用頻域空間域分解法fsdd法分析其模態(tài)參數(shù)。

步驟一中對時(shí)域信號進(jìn)行時(shí)頻分析，如圖2所示，畫出瀑布圖的具體方法如下：

對時(shí)域振動(dòng)信號進(jìn)行stft，得到時(shí)頻譜圖，從所述時(shí)頻譜中可以觀察到階次分量，便于后面步驟的階次提取。

步驟二中由瞬時(shí)頻率估計(jì)法估計(jì)轉(zhuǎn)軸轉(zhuǎn)速的具體方法如下：

如圖3、圖4所示，對振動(dòng)信號進(jìn)行短時(shí)傅里葉變換(short-timefouriertransform,stft)得到時(shí)頻譜，由于最高譜峰能量密度所對應(yīng)的頻率最可能是瞬時(shí)旋轉(zhuǎn)頻率，所以用峰值搜索法對時(shí)頻譜中譜峰能量密度的最大值進(jìn)行提取，從而獲得轉(zhuǎn)軸的旋轉(zhuǎn)頻率，得到轉(zhuǎn)速信號

峰值搜索的過程為：

(1)旋轉(zhuǎn)機(jī)械升降速過程總有一個(gè)穩(wěn)定轉(zhuǎn)速階段，確定穩(wěn)定轉(zhuǎn)速對應(yīng)的頻率值；

(2)以穩(wěn)定時(shí)刻的頻率值為峰值搜索起始點(diǎn)，順次進(jìn)行峰值搜索；

(3)在搜索過程中保證前一時(shí)刻的頻率不大于(升速階段)或不小于(降速階段)后一時(shí)刻的頻率值；

(4)搜索過程中將不滿足條件(3)的特殊點(diǎn)剔除，結(jié)合樣條插值法，這樣就可以得到轉(zhuǎn)軸旋轉(zhuǎn)頻率，從而得到轉(zhuǎn)速信號。

步驟三中自適應(yīng)vold-kalman濾波階比跟蹤算法對某一階次信號進(jìn)行階次提取的具體過程為：

(i)狀態(tài)方程

a(nδt)-2cos(ωδt)a((n-1)δt)+a((n-2)δt)＝0(1)

式中，δt為離散時(shí)間；a(nδt)是第n次離散時(shí)間采樣；ω是正弦波瞬時(shí)頻率。

公式(1)描述了一個(gè)在連續(xù)三個(gè)時(shí)間點(diǎn)頻率幅值都是恒定的正弦波。由于階比的頻率是隨時(shí)間變化的，不是恒定的。狀態(tài)方程可用(2)表示：

a(nδt)-2cos(ωδt)a((n-1)δt)+a((n-2)δt)＝ε(n)(2)

式中：ε(n)稱為非一致項(xiàng)，用來描述理想正弦波幅值和頻率的變化。a(nδt)表示為第n個(gè)采樣點(diǎn)的狀態(tài)，整個(gè)系統(tǒng)的狀態(tài)方程里有n個(gè)采樣點(diǎn)，(2)式展開為：

由式(3)可得狀態(tài)方程得矩陣形式：

fa＝ε(4)

(ii)觀測方程

實(shí)際測得的振動(dòng)信號是由各階比成分的和再加上測量誤差和噪聲組成的，觀測方程描述了階比x(n)和測量數(shù)據(jù)y(n)之間的關(guān)系。測量數(shù)據(jù)不僅包含感興趣的階比，還包含機(jī)器產(chǎn)生的所有階比和背景噪聲。觀測方程表示如下：

y(n)＝x(n)+σ(n)(5)

式中，σ(n)是非跟蹤的階比和隨機(jī)噪聲；x(n)為所提取的階比成分；y((n)為實(shí)測振動(dòng)信號的第n個(gè)采樣點(diǎn)的值。

觀測方程展開：

寫成矩陣形式為：

y＝x+δ(7)

(iii)自適應(yīng)狀態(tài)參數(shù)識別遞推算法

首先要引入一個(gè)離散控制過程的系統(tǒng)，該系統(tǒng)可用一個(gè)線性隨機(jī)微分方程來描述：

a(k)＝fa(k-1)+bu(k)+ψ(k)(8)

再加上系統(tǒng)的測量值：

y(k)＝ca(k)+ξ(k)(9)

上兩式子中，a(k)是k時(shí)刻的系統(tǒng)狀態(tài)，u(k)是k時(shí)刻對系統(tǒng)的控制量。f和b是系統(tǒng)參數(shù)，對于多模型系統(tǒng)，它們?yōu)榫仃?。y(k)是k時(shí)刻的測量值，c是觀測系統(tǒng)觀測矩陣。ψ(k)和ξ(k)分別表示狀態(tài)方程和觀測方程的噪聲。它們被假設(shè)成高斯白噪聲，其協(xié)方差分別是q、r，這里假設(shè)它們不隨系統(tǒng)狀態(tài)變化而變化。

由于滿足上面的條件(線性隨機(jī)微分系統(tǒng)，狀態(tài)和觀測都是高斯白噪聲)，卡爾曼濾波器是最優(yōu)的信息處理器。如圖5所示，下面來估算系統(tǒng)的最優(yōu)化輸出。

首先利用系統(tǒng)的狀態(tài)轉(zhuǎn)移矩陣預(yù)測下一個(gè)狀態(tài)的系統(tǒng)。假設(shè)現(xiàn)在的系統(tǒng)狀態(tài)是k，根據(jù)系統(tǒng)的模型，可以基于系統(tǒng)的上一狀態(tài)而預(yù)測出現(xiàn)在狀態(tài)：

a(k|k-1)＝fa(k-1|k-1)+bu(k)(10)

式(10)中，a(k|k-1)是利用上一個(gè)狀態(tài)預(yù)測的結(jié)果，a(k-1|k-1)是上一個(gè)狀態(tài)最優(yōu)的結(jié)果，u(k)為現(xiàn)在狀態(tài)的控制量，如果沒有控制量，它可以為0。

到現(xiàn)在為止，系統(tǒng)結(jié)果已經(jīng)更新了，可是對應(yīng)于a(k|k-1)的協(xié)方差還沒有更新。用p表示協(xié)方差：

p(k|k-1)＝fp(k-1|k-1)f′+q(11)

式子(11)中p(k|k-1)是a(k|k-1)對應(yīng)的協(xié)方差，p(k-1|k-1)是a(k-1|k-1)對應(yīng)的協(xié)方差，f′表示f的轉(zhuǎn)置矩陣，q是系統(tǒng)狀態(tài)的協(xié)方差。

有了現(xiàn)在狀態(tài)的預(yù)測結(jié)果，再收集現(xiàn)在狀態(tài)的測量值。結(jié)合預(yù)測值和測量值，可以得到現(xiàn)在狀態(tài)(k)的最優(yōu)化估算值a(k|k)：

a(k|k)＝a(k|k-1)+kg(k)(y(k)-ca(k|k-1))(12)

其中kg為卡爾曼增益：

此外，為了卡爾曼濾波器運(yùn)行至系統(tǒng)過程結(jié)束，需更新k狀態(tài)下a(k|k)的協(xié)方差：

p(k|k)＝(i-kg(k)c)p(k|k-1)(14)

其中i為單位矩陣，對于單模型單測量，i＝1。當(dāng)系統(tǒng)進(jìn)入式(k+1)狀態(tài)時(shí)，p(k|k)就是式(11)中的p(k-1|k-1)。

步驟四中對提取的階次信號使用fsdd法分析其模態(tài)參數(shù)具體方法如下：

如圖6所示，提取階次分量后，應(yīng)用fsdd法對階次分量信號進(jìn)行分析識別，具體過程為：求出階次分量的功率譜后，經(jīng)過奇異值分解后得到的左奇異向量與右奇異向量，可以計(jì)算得到一個(gè)增強(qiáng)功率譜，增強(qiáng)功率譜近似于一個(gè)單自由度系統(tǒng)，由此可以計(jì)算得到自然頻率與阻尼比。

奇異值分解公式為：

[u][s][v]^h＝svd({o1}{o2}…{om})(15)

式(15)中，[u]是左奇異向量矩陣；[s]是奇異值對角矩陣；[v]是右奇異向量矩陣；{oi}是階次分量i的功率譜。

基于階次分量的增強(qiáng)功率譜計(jì)算公式為：

g(jω)＝{ul}^h[{o1}{o2}…{om}]n×m{vl}(16)

式(16)中，g(jω)為基于階次分量的增強(qiáng)功率譜；{ul}^h是左奇異向量的共軛轉(zhuǎn)置；{o}是階次分量i的功率譜；{vl}是右奇異向量。