非零奇異值數量的理論分析及其在滑動軸承-轉子振動特征提取應用

2019-08-19 02:08:38楊期江趙學智湯雅連李偉光滕憲斌郭明軍

振動與沖擊 2019年15期

關鍵詞：信號

楊期江，趙學智，湯雅連，李偉光，滕憲斌，郭明軍

(1. 廣州航海學院輪機工程學院，廣州 510725; 2．華南理工大學機械與汽車工程學院，廣州 510640；3. 廣東金融學院互聯網金融與信息工程學院, 廣州 510521)

奇異值分解(Singular Value Decomposition, SVD)是一種正交變換，其理論本質是矩陣的低秩化，使得有用信息集中于非零奇異值序列前部，噪聲信息分布于非零奇異值序列后面，并取非零奇異值序列進行矩陣重構，濾掉噪聲信息，從而達到降噪目的。由于奇異值分解具有以上優點，被廣泛應用于小波變換結果的后續處理[1]，信號消噪[2]、滾動軸承故障診斷[3-4]和早期故障微弱信號分離[5]等方面。

與其他降噪方法相比，SVD具有自身獨特的優勢，其處理結果不存在相位漂移，波形失真小[6]。SVD方法不像小波分析那樣依賴于小波基函數的選擇，也不像經驗模態分解存在端點效應問題及模態混疊現象。但是在應用這一方法時，存在一個關鍵問題需要解決。那就是非零奇異值的選擇直接決定著消噪效果的好壞[7]。如果非零奇異值數目選擇過多，則會使得處理結果混入更多的噪聲，如果非零奇異值數目選擇過少，卻又會丟失信號中的有用頻率成分，造成波形畸變，無法準確恢復我們所需要提取的特定頻率成分的信號。因此針對非零奇異值選擇這一關鍵問題，學者們相繼提出了奇異熵[8]、差分譜理論[9]等，但這些方法并沒用解釋清楚非零奇異值個數背后所蘊含著的真正的物理含義。以振動信號為例，當含有多個頻率成分的確定信號混有噪聲時，采用SVD方法進行消噪，怎么確定非零奇異值的個數，既可以完全恢復原來含有多個頻率成分的確定信號，又能夠準確地提取某個或某幾個頻率成分的信號。錢征文等[10]通過實驗說明了利用信號所構造的Hankel矩陣的秩與信號中的主頻個數存在兩倍的關系，但遺憾的是沒有直接給出這一結果的理論推導及證明過程。趙學智等[11-12]提出了Hankel 矩陣方式下確定性信號的非零奇異值和信號所含頻率數量之間的關系，但沒有從理論上證明非零奇異值和頻率個數之間的規律。

本文通過反復仿真實驗的方法分析了Hankel矩陣下非零奇異值數目與信號中的頻率個數成兩倍的數量關系，驗證了奇異值成對出現規律，當構造的m×n的Hankel矩陣行數與列數充分接近時，信號中同一頻率下的兩個非零奇異值會是緊密排列在一起的。根據Hankel矩陣的構造方式，從理論上證明了非零奇異值數目與單個頻率成2倍的關系，當時間域信號是由多個不同頻率疊加時，非零奇異值數目始終是與頻率個數成2倍的數量關系，且非零奇異值個數與幅值和相位無關。并將這一規律應用于旋轉機械滑動軸承-轉子振動的特征頻率信號提取，實現了對轉子不對中故障軸心軌跡的準確提純。

1 Hankel矩陣下奇異值分解信號消噪理論

當應用SVD進行信號處理時，由于信號是一維的，而SVD針對的是矩陣，因此需要將一維信號向量轉化成矩陣形式才能利用SVD將信號分解到不用的向量空間。因此對一個確定性信號y(t),為了對其進行奇異值分解，首先需要利用y(t)構造一個矩陣，這個矩陣一般為Hankel矩陣，如下式所示：

(1)

式中：y[i]為第i個采樣點，i=0,1,…，N-1,N是y(t)長度，H∈Rm×n，矩陣行數m，列數n和信號長度N之間的約束關系：m=N-n+1。對H進行奇異值分解

H=USVT

(2)

式中:U、V是單位正交矩陣U∈Rm×m及V∈Rn×n,S=[diag(σ1,σ2,…,σq),0]或其轉置，取決于n>m或n

為了說明奇異值降噪的原理，我們將U，V用列向量來表示，U=(u1,u2,…,um)，V=(v1,v2,…,vn)，則式(2)可以改寫為列向量ui和vi的表現形式。

(3)

式中:ui∈Rm×1，vi∈Rn×1，i=1,2,…q,q=min(m,n)，式中共有q個奇異值，但這些奇異值并不都是有效的，對于確定性的無噪聲信號x(t)，它的非零奇異值才是有效奇異值，但是非零奇異值數量不一定就等于q，因為如果沒有噪聲，這q個奇異值中將會有零奇異值存在；而對于含噪聲的信號，在其q個奇異值中，有很多非零奇異值是由噪聲產生的，因此其有效奇異值的數目也不等于q。對于一個具有確定頻率個數的信號，不管它含噪與否，它的非零奇異值數目是與確定的頻率個數相關的。近來有學者針對這一問題進行了實驗上的探索，趙學智等首先發現了非零奇異值數量與信號中的頻率個數成兩倍這一關系，但沒有從理論上證明這一結果。作者經過研究，發現可以從理論上對奇異值的數量規律進行證明。下面先從實驗角度分析奇異值的數量規律，然后給出理論證明過程。

2 Hankel矩陣下非零奇異值數量規律的實驗分析

先通過實驗分析非零奇異值數量與頻率個數之間的關系，采用反復實驗的方法。具體為：從含有一個頻率的信號出發，逐步增加信號中頻率的個數，利用信號構造相同維數的 Hankel 矩陣進行 SVD 處理，觀察并分析非零奇異值數目的變化。采用 8個信號進行試驗，下一個信號在前一個信號的基礎上增加一個頻率，且采樣頻率相同，隨著頻率值逐漸增大，其幅值也逐漸加大，而相位隨機，實驗的8個信號如式(4)所示：

y1(t)=2.8sin(t),

y2(t)=2.8sin(t)+4.6sin(4t+0.56),

y3(t)=2.8sin(t)+4.6sin(4t+0.56)+

6.6sin(12t+1.35)

y4(t)=2.8sin(t)+4.6sin(4t+0.56)+

6.6sin(12t+1.35)+8.4sin(20t+2.58)

y5(t)=2.8sin(t)+4.6sin(4t+0.56)+

6.6sin(12t+1.35)+8.4sin(20t+2.58)+

10.6sin(30t+3.12)

y6(t)=2.8sin(t)+4.6sin(4t+0.56)+

6.6sin(12t+1.35)+8.4sin(20t+2.58)+

10.6sin(30t+3.12)+14sin(60t+2.79)

y7(t)=2.8sin(t)+4.6sin(4t+0.56)+

6.6sin(12t+1.35)+8.4sin(20t+2.58)+

10.6sin(30t+3.12)+14sin(60t+2.79)+

17sin(86t+1.29)

y8(t)=2.8sin(t)+4.6sin(4t+0.56)+

6.6sin(12t+1.35)+8.4sin(20t+2.58)+

10.6sin(30t+3.12)+14sin(60t+2.79)+

17sin(86t+1.29)+21.4sin(104t+3.12)

(4)

設每個信號在[0, 2π]內均采集1 024點，利用采樣點構造512×513的大小的Hankel矩陣，并計算其奇異值，得到各個信號的奇異值結果具體如圖1所示，圖中從左至右依次為y1(t)～y8(t)的奇異值曲線。從圖 1可見，當信號每增加一個頻率成分，不管這個幅值、相位如何，非零奇異值的數目總是只增加兩個。在下一個信號的奇異值曲線中，除了新增加的兩個奇異值較大外，其余的奇異值均和上一個信號的奇異值大小幾乎相等，而且新增的兩個奇異值之間幾乎相等。采用反復實驗方法，在仿真信號(4)中，如果采用其它的頻率和幅值，也會得到同樣結果：即不管頻率大小及其幅值如何，一個頻率總是只產生兩個非零奇異值。

為了分析新增的兩個奇異值之間幾乎兩兩相等的原因。我們構造一個含有4個頻率的信號，具體如式(5)所示：

y(t)=sin(5t)+2.6sin(9t+0.35)+

3.7sin(19t+1.26)+5.4sin(32t+0.98)

(5)

對此信號在[0, 2π]內均采集1 024點，分別構造不同維數大小的Hankel 矩陣，并計算其奇異值。在不同維數的Hankel矩陣下各個信號的奇異值結果具體如圖2所示。從圖 2可見，隨著構造的m×n的Hankel矩陣維數增加，奇異值大小都有不同程度的增大。當行數與列數充分接近時，會出現信號中同一頻率下的兩個奇異值幾乎相等。

圖2 非零奇異值與矩陣維數之間的關系

通過上述方法的反復實驗，我們發現了非零奇異值數目與頻率個數之間的數量關系規律以及奇異值成對出現規律：① 對于一個含有固定頻率數量的確定性信號，利用其構造m×n的 Hankel 矩陣，當m和n中最小值大于信號中頻率數量的兩倍之后，非零奇異值的數目始終為信號中頻率數量的兩倍，稱之為奇異值數量規律。② 當構造的m×n的Hankel矩陣行數與列數充分接近時，會出現信號中同一頻率下的兩個非零奇異值是緊密排列在一起的，即非零奇異值成對出現規律。Hankel矩陣下非零奇異值與信號中頻率成分個數之間呈2倍的奇妙的數量關系規律。通過研究，我們發現這一現象背后存在理論上的必然性，下面從理論上來證明這一規律。

3 非零奇異值數量規律的理論證明

設有時間域連續的確定性信號y(t)包含I個頻率成分，由下式表示：

(6)

以采樣間隔為Ts，采樣點數為j，那么時間離散為t=Tsj，因此將原始連續信號離散為y[j]。根據y[j]構建一個m×n的Hankel矩陣H，顯然H可以表示為單個頻率分量構造的Hankel矩陣Hi之和[1]，即：

(7)

式中Hi是頻率分量yi[j]=aisin(ωiTs·j+φi)構造的m×n的Hankel矩陣，如下式所示[2]：

(8)

式中:m=N-n+1，而N是xi(t)離散后的長度。

現在分三種情況分別對Hankel矩陣H進行討論，第一種情況為確定性信號y(t)僅含有1個頻率成分，第二種情況為確定性信號y(t)含有I個頻率成分，第三種情況為確定性信號y(t)含有噪聲時。

(1) 當確定性信號y(t)僅含有1個頻率成分時

此時原始信號y(t)僅表示為：

y(t)=A1sin(ω1t+φ1)

(9)

那么其離散之后構建的m×n的Hankel矩陣為：

(10)

根據三角函數關系式，可得：

sin(ω1t+φ1)=a1sin(ω1t)+b1cos(ω1t)

(11)

將式(11)代入式(10)，可得：

(12)

由歐拉公式ejω1Ts=cosω1Ts+jsinω1Ts，得：

(13)

將式(13)代入式(12)中得：

(14)

通過化簡為：

(15)

可將H1改寫為兩個Hankel矩陣之和的形式：

(16)

其中：

(2) 當確定性信號y(t)含有I個頻率成分時

此時原始信號y(t)表示含單個頻率時域信號之和，即為：

(17)

那么其離散之后構建的m×n的Hankel矩陣H為單個頻率分量構造的Hankel矩陣Hi之和：

(18)

同理，將Hi改寫為兩個Hankel矩陣之和的形式：

(19)

式中：

將式(19)代入到式(18)中，可得含有I個頻率成分信號總Hankel矩陣可寫成如下形式：

(20)

由式(20)可知，矩陣H由I個Hankel矩陣Hi求和得到，其中第i個頻率成分信號的Hankel矩陣Hi的秩為2，僅與頻率ωi有關，與幅值Ai，相位角φi都無關。由于ω1≠ω2≠ω3≠…≠ωI，可知每一個頻率信號所構成的Hankel矩陣Hi之間是互相線性無關的，因此可知當原始信號含有I個頻率成分時，總矩陣H的秩為2I，當且僅當m行×n列的總矩陣H滿足min(m,n)≥2I時，總矩陣H秩才為2I。由非零奇異值數量等于矩陣的秩，因此I個頻率成分信號構造的m行×n列Hankel矩陣滿足min(m,n)≥2I時，通過SVD分解之后的非零奇異值個數也應該為2I個，且非零奇異值個數與幅值和相位無關。

(3) 當確定性信號y(t)中含有噪聲時

噪聲的頻譜是分布在整個頻率軸上的，因此噪聲將產生很多非零奇異值，但是噪聲產生的非零奇異值只是影響確定性信號的非零奇異值大小，而不會影響確定性信號的非零奇異值個數。而且非零奇異值的大小是由頻率的幅值決定的，某頻率的幅值越大，其對應的非零奇異值也越大，因為奇異值序列是降序排列的，如果確定性信號的幅值較大，經過SVD處理后，噪聲產生的非零奇異值就會排列在確定性信號的非零奇異值后面，因此，根據確定性信號中的頻率數量，選擇前面合適的奇異值進行SVD重構，就可以消除噪聲的影響，得到干凈的信號。

綜上理論分析可得到如下重要結論：采用Hankel矩陣構造一維信號，如果一維信號中所含互不相等的頻率成分個數為I，當矩陣當中的行數與列數的最小值大于頻率成分個數的兩倍時，通過SVD之后的非零奇異值個數應該為2I個，可表述為每兩個非奇異值對應信號中的一個頻率成分；且非零奇異值個數是與幅值和相位無關的，即非零奇異值數量規律。

對一維信號所構造的Hankle矩陣進行奇異值分解，得到從大到小依次排列的非零奇異值序列譜。根據以上的實驗及理論分析可知，振動信號的特征頻率成分就隱含在這一系列的非零奇異值序列譜當中。下面通過非零奇異值和頻率的關系來提取滑動軸承-轉子振動信號特征頻率。

4 對滑動軸承-轉子振動信號特征的提取

滑動軸承是大重型旋轉機械的核心支撐部件。軸承間隙中所形成的楔形油膜支撐轉子高速旋轉。相對于滾動軸承，滑動軸承與轉子之間存在便于形成油膜的間隙，大型機組在進行安裝調試時，聯軸器與轉子之間的對中精度往往是制約轉子高速運行的關鍵因素。因此機組在進行現場安裝前，需要利用實驗臺調試模擬實際安裝環境，通過對滑動軸承-轉子進行振動測試提取軸心軌跡等特征信號以判斷滑動軸承-轉子的運行狀況。實際測試過程中，機組試驗臺的振動受到各個部件耦合振動的影響，所測試得到的實驗數據噪聲干擾較大，導致軸心軌跡等特征信號模糊不清，難以辨別。

(a) 試驗臺架實物圖

(b) 電渦流傳感器安裝圖

本文試驗研究采用的軸承-轉子試驗臺為團隊自主研發搭建的大型滑動軸承轉子試驗臺系統。最高轉速5 000 r/min；轉子重量約5T，具體如圖3(a)所示。圖3(b)是實際機組試驗臺電渦流位移傳感器安裝圖。本文轉子試驗轉速為2 600 r/min，圖4是滑動軸承-轉子近電機端實際運行過程中采集到的軸心軌跡，通過該軸心軌跡難以判斷機組試驗臺運行狀態。

圖4 滑動軸承-轉子振動軸心軌跡

為了更進一步地分析影響軸心軌跡的各個頻率成分，首先分別對近電機端1號測點(x方向)與2號測點(y方向)上采集得到位移振動原始信號進行快速傅里葉變化(FFT)，可得到其頻率幅值譜如圖5所示。從圖5(b)，(d)頻域圖可知，1倍頻與2倍頻幅值最為突出，且2倍頻與1倍頻幅值接近，同時存在高倍轉頻，根據轉子振動理論，二倍轉頻的存在并伴有高倍轉頻這種現象主要是由于轉子對中不良所致，其中還有明顯的電力系統的工頻干擾，分別是50 Hz與150 Hz。如果只通過單個測點信號的頻譜分析難以判別轉子不對中的程度，需要引入軸心軌跡以進一步識別，而由于高倍轉頻以及電力工頻等頻率成分疊加導致原始信號所繪制的軸心軌跡模糊，難以辨識。

根據Hankel矩陣下非零奇異值數目與頻率個數之間的數量關系，以及奇異值成對出現規律，對原始信號1倍頻與2倍頻成分進行特征提取。原信號的數據長度為4 096，故利用圖5(a)、(c)的時域信號構造兩個2 048×2 049的Hankel矩陣并對其進行奇異值分解，得到兩個原信號的奇異值如圖6所示(圖中只列出前40個奇異值)。

原信號中主要有1倍頻與2倍頻這兩個幅值較大的頻率成分，它們產生4個奇異值，而由于頻率的幅值對奇異值的影響最大[11]，因此這2個幅值較大的頻率對應的很可能就是前4個奇異值，現在利用第1對和第2對奇異值進行 SVD 重構，即有重構矩陣為[12]：

(21)

(a) 近電機端1號測點時域圖

(b) 近電機端1號測點頻域圖

(d) 近電機端2號測點頻域圖

圖5 滑動軸承轉子振動信號及其幅值譜

Fig.5 Vibration signal of sliding bearing-rotor and its spectrum

(a) 近電機端1號測點

(b) 近電機端2號測點

得到的矩陣H′也是一個2 048×2 049的矩陣，采用平均法[12]從H′中恢復出信號，結果如圖7(b)、(d)所示。從幅值譜可見這兩個重構信號的頻率成分確實為轉頻和二倍轉頻，它們的幅值分別為1.542×10-5m和1.586×10-5m、1.266×10-5m和2.152×10-5m，而在兩個原信號的幅值譜中 (圖5(b)、(d)所示)，這兩個頻率的幅值分別為1.539×10-5m和1.594×10-5m、1.269×10-5m和2.15×10-5m，SVD提取的結果與它們只存在微小的差異，這種差異是噪聲的影響造成的，是不可避免的，可以說，SVD把轉頻和二倍轉頻基本完整地提取出來了。從提取出的時域波形(圖7(a)、(c)所示)可以看出1倍頻與二倍頻所對應的時域波形都沒有發生較大的波動，因此轉子不對中狀態是比較穩定的。從提取出的頻域波形(圖7(b)、(d)所示)可以看出該頻域圖是非常干凈的，只有1倍轉頻和2倍頻兩個譜峰存在，尤其令人注意的是，不但這兩個譜峰之間的電力系統50 Hz的工頻干擾也完全消失了，而且達到了對高次諧波干擾完美的濾波效果，這種效果是一般帶通濾波器達不到的?？梢钥闯鯯VD分離特定頻率成分的效果是非常優異的，一般的濾波器很難做到這一點，濾波器總是會存在一定的帶寬，并且還有過渡帶的影響，除非兩個頻率隔得特別遠，否則濾波器在分離單個頻率時總是會受到臨近頻率的影響。

(a) 近電機端1號測點時域圖

(b) 近電機端1號測點頻域圖

(d) 近電機端2號測點頻域圖

圖7 SVD提取的故障波形特征及其幅值譜

Fig.7 Waveform feature of fault extracted by SVD and its spectrum

分別以圖7(a)、(c)時域波形為x軸、y軸得到滑動軸承-轉子近電機端振動軸心軌跡，具體如圖8(a)所示。同理，采用上述方法可以得到滑動軸承-轉子遠電機端振動軸心軌跡如圖8(b)所示。根據轉子不對中故障的特征分析，當不對中比較輕微時，軸心軌跡呈橢圓狀；當不對中故障達到中等程度時，軸心軌跡呈香蕉形；當不對中故障較嚴重時，軸心軌跡呈外“8”字形。從圖8可以看出，滑動軸承-轉子近電機端振動軸心軌跡呈香蕉狀，不對中故障為中等程度；滑動軸承-轉子遠電機端振動軸心軌跡呈外“8”字形，不對中故障較嚴重，所以遠電機端的不對中程度要大于近電機端。

(a) 近電機端

(b) 遠電機端

經驗模態分解(Empirical Mode Decomposition，EMD)是一種廣泛應用的信號分解方法，也可以分離出信號中的各個頻率。近年來發展了EEMD(總體平均經驗模態分解)、CEEMD(補充的 EEMD)等不同的改進方法[13]。對于本文的滑動軸承-轉子振動信號，利用CEEMD方法自適應分解得到一系列從高頻到低頻的具有不同物理意義內稟模態分量(Intrinsic Mode Function, IMF)，具體如圖9(a)、(c)所示。通過快速傅里葉變換(FFT)得到各IMF的幅值譜，具體如圖9(b)、(d)所示。從圖9(b)可知近電機端1號測點信號經CEEMD方法分解之后得到的IMF3的頻譜圖中主要是以2倍轉頻為主，但是沒有能夠分離電力系統的工頻干擾成分，如1倍工頻50 Hz，3倍工頻150 Hz；IMF4主要是以1倍轉頻為主，但是仍然存在模態混疊現象，如2倍頻的存在。從圖9(d)可知近電機端2號測點信號CEEMD分解之后得到的IMF3的頻譜圖中主要是以2倍轉頻為主，但同樣還伴有電力系統的工頻干擾，如1倍工頻50 Hz，3倍工頻150 Hz，同時仍然存在模態混疊現象，如1倍頻的存在；IMF4主要是以1倍轉頻為主。因此利用內稟模態分量IMF3與IMF4進行CEEMD重構，即有重構信號為：

(22)

式中：ci為第i個IMF。

可得到近電機端1號、2號測點提取的故障波形，分別以這兩個重構之后的時域波形為x軸、y軸得到滑動軸承-轉子近電機端振動軸心軌跡，具體如圖10(a)所示。同理，采用上述方法可以得到滑動軸承-轉子遠電機端振動軸心軌跡如圖10(b)所示。與圖8對比可見，SVD方法提取的軸心軌跡比EMD方法提取的軸心軌跡清楚很多。