基于M-K檢驗(yàn)和最小二乘法的SSA算法研究①

賈國(guó)慶, 程明龍,*, 易輝躍, 房衛(wèi)東, 許 暉

(1.青海民族大學(xué)物理與電子信息工程學(xué)院，青海 西寧 810007；2.中國(guó)科學(xué)院上海微系統(tǒng)與信息技術(shù)研究所，上海 201899；3.上海無(wú)線通信研究中心，上海 201210)

0 引 言

通過(guò)振動(dòng)信號(hào)獲取結(jié)構(gòu)體的固有頻率是監(jiān)測(cè)結(jié)構(gòu)體的工作狀態(tài)及健康狀況常用的方式。在振動(dòng)信號(hào)的采集過(guò)程中，會(huì)夾雜一些直流分量和噪聲，對(duì)分析結(jié)果產(chǎn)生影響；同時(shí)，傳感器溫度造成的零點(diǎn)漂移會(huì)產(chǎn)生線性趨勢(shì)項(xiàng)，使信號(hào)嚴(yán)重失真[1-2]。如何消除不利影響因素一直是振動(dòng)信號(hào)處理中的研究重點(diǎn)。文獻(xiàn)[3]結(jié)合自適應(yīng)噪聲完全集合經(jīng)驗(yàn)?zāi)B(tài)分解和奇異譜分析，對(duì)含噪信號(hào)進(jìn)行自適應(yīng)分解，然后使用SSA對(duì)信號(hào)層進(jìn)一步降低低頻噪聲，能夠較好地消除噪聲的干擾。文獻(xiàn)[4]提出一種基于相位空間重構(gòu)和奇異譜分析的降噪法，通過(guò)奇異值分解對(duì)比噪聲信號(hào)和有用信號(hào)的奇異值數(shù)目進(jìn)行自適應(yīng)降噪，重構(gòu)原始數(shù)據(jù)。文獻(xiàn)[5]介紹了一種滑動(dòng)窗長(zhǎng)SSA算法，通過(guò)改變軌跡矩陣的窗長(zhǎng)，突破了特征分解的局限性，更利于提取真實(shí)信號(hào)。上述算法對(duì)信號(hào)突變嚴(yán)重的情況下有一定的優(yōu)勢(shì)，但是相位移動(dòng)和奇異值突變的現(xiàn)象沒(méi)有得到改進(jìn)。針對(duì)上述算法的不足，提出了優(yōu)化算法，通過(guò)對(duì)信號(hào)添加偽序列，并融合曼-肯德?tīng)?Mann—Kendall, M-K)檢驗(yàn)法和最小二乘法，可以避免信號(hào)相移的出現(xiàn)，降低突變奇異值的能量，自適應(yīng)的選擇分組重構(gòu)階次，從而增加SSA的降噪能力。仿真結(jié)果表明，改進(jìn)SSA算法可以達(dá)到更好的降噪效果，提高了信號(hào)的信噪比，從而有效提取速度、位移和固有頻率等信息。

1 奇異譜分析

奇異譜分析是研究非線性時(shí)間序列的高效方式，通過(guò)構(gòu)建軌跡矩陣、SVD分解、分組和重構(gòu)，可以獲得原時(shí)間序列的趨勢(shì)信號(hào)、周期信號(hào)和噪聲信號(hào)，進(jìn)而對(duì)序列進(jìn)行分析，以去除噪聲信號(hào)。

1.1 軌跡矩陣

存在有限長(zhǎng)度的一維時(shí)間序列[x1,x1,…,xN]，其中N為序列長(zhǎng)度。選擇一個(gè)合適的窗口長(zhǎng)度L,將原始一維時(shí)間序列進(jìn)行滯后排列，得到軌跡矩陣[6]:

(1)

窗口長(zhǎng)度L通常小于N/2。令K=N-L+1，則軌跡矩陣X為：

(2)

1.2 奇異值分解

特征值的分解可以獲取特征值與特征向量，以此分離振動(dòng)信號(hào)的趨勢(shì)項(xiàng)和周期項(xiàng)。在SSA中，使用SVD方式對(duì)軌跡矩陣進(jìn)行奇異值分解式(3)：

X=UΣVT

(3)

式(3)中，U為左奇異值矩陣；Σ為對(duì)角矩陣，即奇異值;V為右奇異值矩陣。U和V均為單位正交矩陣[7]，滿足UUT=I,VVT=I，令

S=XXT

(4)

對(duì)S進(jìn)行奇異值分解，得到對(duì)應(yīng)奇異值λ和奇異值向量U=[U1,U2,…,UL]。對(duì)λ進(jìn)行降序排列：λ1>λ2>…>λr?λr+1>…λL≥0，第r個(gè)奇異值遠(yuǎn)大于第r+1個(gè)奇異值，說(shuō)明在第r個(gè)奇異值之后的奇異值即為噪聲信號(hào)?？梢愿鶕?jù)重新排序后的奇異值的能量占比情況，選擇含有有效信號(hào)的前r個(gè)奇異值向量。

對(duì)于實(shí)際信號(hào)序列有式(5)：

(5)

式(5)中，為確保窗口長(zhǎng)度小于N/2，d=min{L,K}。所以,軌跡矩陣X的SVD分解結(jié)果為式(6)：

X=X1+X2+…+Xd

(6)

1.3 分組和重構(gòu)

將式(6)的下標(biāo)d劃分為m個(gè)互不相交的子集I1,I2,…,Im，令I(lǐng)={i1,i2,…,ip}，則對(duì)應(yīng)I的矩陣XI=Xi1+Xi2+…+Xip。因此可得式(7)：

X=XI1+XI2+…+XIm,(1≤m≤L)

(7)

為了得到分解后的序列，這里將式(7)中的每個(gè)矩陣XIj轉(zhuǎn)換為長(zhǎng)度為N的序列。設(shè)定Y為L(zhǎng)×K維的矩陣，元素為Aij,1≤i≤L,1≤j≤K。令Lp=min(L,K),Kp=max(L,K)，若L

(8)

2 SSA算法優(yōu)化

2.1 數(shù)據(jù)重組

在進(jìn)行傳統(tǒng)SSA算法測(cè)試的過(guò)程中發(fā)現(xiàn)，重組后的信號(hào)在起始和結(jié)束兩端具有相位移動(dòng)的現(xiàn)象，會(huì)對(duì)速度信號(hào)和位移信號(hào)的轉(zhuǎn)換造成一定的影響，使信號(hào)出現(xiàn)失真現(xiàn)象。在不同信噪比和窗口長(zhǎng)度的情況下進(jìn)行測(cè)試得出：信號(hào)兩端相移端點(diǎn)的長(zhǎng)度M約為SSA窗長(zhǎng)L的1/10?；诖?，可以構(gòu)建兩段長(zhǎng)度為M的偽數(shù)據(jù)xp1和xp2，分別添加在待測(cè)數(shù)據(jù)的兩端，使得原始數(shù)據(jù)長(zhǎng)度變?yōu)長(zhǎng)+2M，然后再進(jìn)行重組數(shù)據(jù)的SSA分解重構(gòu)。

X+=[xp1,X,xp2]

(9)

式(9)中，xp1和xp2可以為同一段數(shù)據(jù)。為了使偽數(shù)據(jù)更切合原始數(shù)據(jù)的特征，復(fù)制原始數(shù)據(jù)的前M個(gè)數(shù)據(jù)作為偽數(shù)據(jù)進(jìn)行重組。

2.2 奇異值擬合

當(dāng)信號(hào)的信噪比較低且信號(hào)存在高強(qiáng)度尖峰噪聲時(shí)，SVD并不能完全分離噪聲信號(hào)和有效信號(hào)，通常奇異值矩陣中會(huì)含有突變的奇異值，即噪聲特征被當(dāng)作趨勢(shì)特征。

為解決奇異值突變問(wèn)題，對(duì)分解后的奇異值矩陣進(jìn)行曲線擬合，進(jìn)一步消除殘余噪聲。有式(10)：

(10)

在曲線擬合時(shí)，按偏差平方和最小原則選取擬合曲線式(11)：

(11)

若使δi值最小，需依次求δi對(duì)aj的偏導(dǎo)式(12)：

(12)

2.3 分組重構(gòu)階次的選擇

在對(duì)奇異值降序排列后發(fā)現(xiàn)，通過(guò)選取前r個(gè)能量占比較高的奇異值進(jìn)行信號(hào)重構(gòu)，會(huì)取得較好的降噪效果，r值的選擇通常根據(jù)主觀經(jīng)驗(yàn)進(jìn)行選取，而r值的大小會(huì)對(duì)信號(hào)的處理結(jié)果造成很大的影響。改進(jìn)的SSA算法中，使用M-K檢驗(yàn)法自適應(yīng)選擇r值。

對(duì)X+進(jìn)行SVD分解后可得奇異值向量U，令Uj=(U1,U2,…,UL-1)，有式(13)：

(13)

式(13)中sgn為檢驗(yàn)函數(shù)，令θ=Xk-Xj，如式(14)：

(14)

然后計(jì)算S的方差和M-K檢驗(yàn)統(tǒng)計(jì)量Z[8]：

孔徑分布和截留相對(duì)分子量采用一系列的PEG截留率結(jié)果并通過(guò)Matlab軟件采用下述對(duì)數(shù)正態(tài)分布方程模擬計(jì)算得到[15]：

(15)

(16)

根據(jù)檢驗(yàn)統(tǒng)計(jì)值可以得知：當(dāng)Z>0時(shí)，檢驗(yàn)序列呈上升趨勢(shì);當(dāng)Z<0時(shí),檢驗(yàn)序列呈下降趨勢(shì)；當(dāng)Z=0時(shí)，檢驗(yàn)序列保持原始數(shù)據(jù)的趨勢(shì)，則可以根據(jù)此時(shí)r值的大小進(jìn)行信號(hào)重構(gòu)。

2.4 改進(jìn)算法的流程

SSA改進(jìn)算法的計(jì)算流程如下：

(1)選定窗口長(zhǎng)度L，確定偽數(shù)據(jù)長(zhǎng)度M，分別添加至原始數(shù)據(jù)兩端；

(2)對(duì)重組的數(shù)據(jù)進(jìn)行滯后排列，創(chuàng)建軌跡矩陣；

(3)分解軌跡矩陣，獲得奇異值和奇異值矩陣，并對(duì)奇異值降序排列；

(4)對(duì)奇異值矩陣進(jìn)行曲線擬合，去除突變奇異值；

(5)預(yù)設(shè)定重構(gòu)階次r=1，進(jìn)行趨勢(shì)檢驗(yàn)，若檢驗(yàn)統(tǒng)計(jì)量Z值不為0，則r=r+1，直至Z值為0為止；

(6)將軌跡矩陣映射至包含前r個(gè)奇異值的奇異值矩陣上，計(jì)算時(shí)間主成分；

(7)求時(shí)間主成分矩陣的對(duì)角線平均，獲得重構(gòu)后的振動(dòng)加速度信號(hào)。

3 實(shí)驗(yàn)仿真

利用計(jì)算機(jī)仿真驗(yàn)證所提改進(jìn)SSA算法的有效性。實(shí)驗(yàn)中使用正弦信號(hào)，采樣頻率fs=2000Hz，采樣點(diǎn)數(shù)N=1024。在信噪比為-20dB～20dB的環(huán)境下進(jìn)行30000次的蒙卡特羅仿真模擬。假設(shè)加速度信號(hào)為含有高斯白噪聲的正弦信號(hào)，可表示為式(17)：

s(t)=A×sin(2×π×f×t)+n(t)

(17)

式(17)中，f為信號(hào)頻率；n(t)為均值為0的高斯白噪聲。

進(jìn)行SSA分析之前，先消除因傳感器溫度等因素造成的零點(diǎn)漂移式(18)：

(18)

仿真中，f=20Hz，窗口長(zhǎng)度L=300，偽數(shù)據(jù)長(zhǎng)度M=30。SSA分析結(jié)果如圖1所示，其中SNR=-4dB。對(duì)噪聲信號(hào)依次進(jìn)行了傳統(tǒng)SSA，滑動(dòng)窗長(zhǎng)SSA(SSA-L)，基于EEMD的SSA(SSA-E)和本文所提的改進(jìn)SSA處理。

由圖1可知，由于SSA強(qiáng)大的分解能力，經(jīng)其處理后的信號(hào)不會(huì)存在諧波分量，更加平滑。傳統(tǒng)的SSA算法和SSA-L算法，SSA-E算法所重構(gòu)出的信號(hào)都存在相位移動(dòng)；由于分解過(guò)程中的特征誤選，導(dǎo)致部分信號(hào)能量提高或降低。而改進(jìn)的SSA算法通過(guò)偽數(shù)據(jù)有效地解決了處理后的加速度信號(hào)在起始和結(jié)束兩端的相位移動(dòng)現(xiàn)象，改進(jìn)后的奇異值能量趨于正常，誤選幾率降低，奇異值矩陣能夠最大程度表述加速度信號(hào)的特征，使得結(jié)果更加逼近理想信號(hào)。

為了驗(yàn)證所提改進(jìn)SSA算法的去噪效果，在信噪比為-20dB～20dB的環(huán)境下，求取了降噪后信號(hào)的信噪比如式(19)：

(19)

式(19)中，Ps為加速度信號(hào)的能量；Pn為噪聲信號(hào)的能量。各算法降噪后的信噪比對(duì)比結(jié)果如圖2所示。SSA分解中的奇異值突變，會(huì)導(dǎo)致部分噪聲奇異值能量不能迅速衰減，從而被當(dāng)作加速度特征進(jìn)行分組。改進(jìn)的SSA算法通過(guò)曲線擬合提高了有效信號(hào)的能量，同時(shí)優(yōu)化了奇異值個(gè)數(shù)r的選擇，極大地提高了信號(hào)的質(zhì)量，尤其是在低信噪比的情況下，效果更為顯著。

降噪質(zhì)量的高低決定著是否能夠準(zhǔn)確獲取相關(guān)的速度和位移等信息。對(duì)加速度信號(hào)進(jìn)行一次積分和二次積分分別得到速度信號(hào)和位移信號(hào)。如圖3和圖4分別給出了SNR=-4dB時(shí)由加速度信號(hào)積分得到的速度信號(hào)和位移信號(hào)。由圖3和圖4可見(jiàn)，無(wú)論速度還是位移，改進(jìn)SSA算法更逼近于無(wú)噪聲環(huán)境下的速度和位移信號(hào)。

為了更直觀地表述降噪后加速度信號(hào)的轉(zhuǎn)換效果，使用式(20)求出了速度信號(hào)和位移信號(hào)的均方誤差：

(20)

最后，驗(yàn)證利用所提改進(jìn)SSA算法由位移信號(hào)估計(jì)頻率的正確性。對(duì)位移信號(hào)進(jìn)行FFT，利用頻譜峰值進(jìn)行固有頻率的估計(jì)，估計(jì)值如表1所示。在改進(jìn)的SSA算法下，由于極大地提高了信號(hào)的信噪比，使得頻率估計(jì)值和信號(hào)原始振動(dòng)頻率一致，低信噪比環(huán)境中的估計(jì)結(jié)果優(yōu)于其他的SSA算法。說(shuō)明使用改進(jìn)的SSA算法降噪后的振動(dòng)信號(hào)，攜帶的噪聲分量更低，更利于估計(jì)結(jié)構(gòu)體的固有頻率，從而更有效地對(duì)結(jié)構(gòu)體的工作狀態(tài)和健康狀況進(jìn)行分析。

表1 頻率估計(jì)(Hz)

4 結(jié) 語(yǔ)

傳統(tǒng)的SSA算法在對(duì)原始加速度信號(hào)分解重構(gòu)后，會(huì)存在相位移動(dòng)現(xiàn)象，造成數(shù)據(jù)起始和結(jié)束兩端失真，進(jìn)而影響速度和位移信號(hào)的準(zhǔn)確估計(jì)，并且在關(guān)鍵的重構(gòu)階次選擇上沒(méi)有選值標(biāo)準(zhǔn)，從而對(duì)信號(hào)的降噪效果產(chǎn)生較大的影響?；贛-K檢驗(yàn)和最小二乘法的SSA改進(jìn)算法針對(duì)傳統(tǒng)SSA算法的缺點(diǎn)加以改進(jìn)，同時(shí)保留了傳統(tǒng)算法高效分解信號(hào)的優(yōu)點(diǎn)。在相同的環(huán)境下，利用計(jì)算機(jī)仿真與傳統(tǒng)SSA和其他改進(jìn)算法SSA-L，SSA-E等降噪算法做了對(duì)比。仿真結(jié)果表明，基于M-K檢驗(yàn)和最小二乘法的SSA改進(jìn)算法在任何信噪比的環(huán)境下均優(yōu)于其他幾種算法，對(duì)建筑、橋梁等結(jié)構(gòu)體的振動(dòng)監(jiān)測(cè)診斷提供了有效的信號(hào)處理方式。