黃海陽，王 成+，李海波，賴雄鳴，緱 錦，張惠臻

(1.華僑大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，福建 廈門 361021;2.華僑大學(xué) 機(jī)電及自動(dòng)化學(xué)院，福建 廈門 361021)

0 引言

模態(tài)參數(shù)(模態(tài)固有頻率、模態(tài)振型、模態(tài)阻尼比)能夠很好地反映結(jié)構(gòu)的動(dòng)力學(xué)特性，因此提取模態(tài)參數(shù)成為結(jié)構(gòu)系統(tǒng)辨識(shí)問題的核心[1]。然而，許多大規(guī)模的工程結(jié)構(gòu)輸入無法測(cè)量或代價(jià)昂貴，此時(shí)工作模態(tài)分析方法(又稱僅輸出辨識(shí)方法)能夠在工作環(huán)境中僅從輸出響應(yīng)中提取結(jié)構(gòu)的模態(tài)參數(shù)，因此得到廣泛關(guān)注[2]。目前，工作模態(tài)分析已經(jīng)廣泛應(yīng)用在故障診斷[3]、結(jié)構(gòu)優(yōu)化設(shè)計(jì)[4-5]等方面。

現(xiàn)實(shí)中工程結(jié)構(gòu)大多是時(shí)變的，物理特性(質(zhì)量、剛度、阻尼等)隨時(shí)間而變化，其模態(tài)參數(shù)也會(huì)隨著時(shí)間而變化[6]，例如，火車過橋，火箭發(fā)射，旋轉(zhuǎn)機(jī)械等都呈現(xiàn)出時(shí)變特性[7]。目前，線性時(shí)變結(jié)構(gòu)的模態(tài)參數(shù)識(shí)別方法主要分為時(shí)域法和頻域法(時(shí)頻分析法)，馬志賽等[8]和Zhou等[9]對(duì)這些方法進(jìn)行了綜述，并總結(jié)了每種方法的優(yōu)缺點(diǎn);振動(dòng)響應(yīng)信號(hào)也不能夠提前一次性獲得，而是隨著時(shí)間慢慢采樣得到，因此需要在線識(shí)別方法[10];滑動(dòng)窗方法是一種在線識(shí)別慢時(shí)變系統(tǒng)模態(tài)參數(shù)的方法，已經(jīng)應(yīng)用在一些算法上[11];滑動(dòng)窗方法相比于遞推的方法更能跟蹤時(shí)變特性，因?yàn)槠涮砑恿诵聰?shù)據(jù)同時(shí)舍棄了舊數(shù)據(jù)，減輕了歷史數(shù)據(jù)的影響[12]。

盲源分離(Blind Source Separation，BSS)是在源信號(hào)和混合方式未知的情況下，僅從觀測(cè)信號(hào)中恢復(fù)源信號(hào)的方法[13]，包括獨(dú)立成分分析、稀疏成分分析、二階盲辨識(shí)等[14]。Sadhu等[15]對(duì)BSS方法在工作模態(tài)參數(shù)識(shí)別上的應(yīng)用進(jìn)行了綜述，并指出幾乎所有基于BSS的模態(tài)識(shí)別方法都是在離線或批處理模式下進(jìn)行的，不適用于控制應(yīng)用和參數(shù)的自動(dòng)化識(shí)別。但識(shí)別時(shí)變結(jié)構(gòu)的模態(tài)參數(shù)，需要在線的自適應(yīng)的BSS方法,Luo等[16]對(duì)無線自適應(yīng)處理BSS進(jìn)行了綜合研究;Cardoso等[17]提出了基于連續(xù)更新思想的等變自適應(yīng)源分離(Equivariant Adaptive Separation via Independence, EASI)方法，該方法采用相對(duì)梯度，具有等變化性，迭代的在線更新每個(gè)時(shí)刻的分離矩陣;Amini等[18]將EASI應(yīng)用到線性時(shí)不變結(jié)構(gòu)的工作模態(tài)在線辨識(shí)中，并指出EASI具有收斂性依賴于源統(tǒng)計(jì)量和分析復(fù)雜的缺點(diǎn)。EASI是一種在線方法，能夠節(jié)省存儲(chǔ)空間，在數(shù)據(jù)輸入的同時(shí)就進(jìn)行處理，但是其使用的固定步長不能兼顧收斂性和跟蹤能力[19]。為此，研究者們提出了大量的變步長EASI算法改進(jìn)EASI算法的收斂性和穩(wěn)定性。陸建濤等[20]提出基于分離指標(biāo)的變步長EASI算法，利用一個(gè)非線性單調(diào)遞增函數(shù)自適應(yīng)調(diào)節(jié)步長;鄔金松等[21]探索了分離程度的概念,并根據(jù)信號(hào)的分離程度自適應(yīng)調(diào)整EASI算法的步長。但這些變步長方法都只是在時(shí)不變系統(tǒng)中提高算法性能，時(shí)變系統(tǒng)需要更復(fù)雜的變步長方法。目前，最新改進(jìn)的變步長EASI算法在系統(tǒng)突變后能夠自適應(yīng)追蹤并且重新收斂于新的平衡點(diǎn)[22-23],但當(dāng)系統(tǒng)參數(shù)和源信號(hào)隨時(shí)間持續(xù)變化時(shí)，變步長EASI算法也難以有效跟蹤其特性。

綜上所述，本文將滑動(dòng)窗方法[24]與變步長EASI算法相結(jié)合，提出一種基于滑動(dòng)窗變步長等變自適應(yīng)源分離(Moving Window Variable step-size EASI, MWVEASI)的線性慢時(shí)變結(jié)構(gòu)工作模態(tài)參數(shù)在線識(shí)別方法。利用滑動(dòng)窗的在線跟蹤特性和變步長EASI方法的收斂性，識(shí)別線性時(shí)變結(jié)構(gòu)的瞬態(tài)模態(tài)固有頻率和振型。

1 相關(guān)工作

1.1 傳統(tǒng)盲源分離問題描述

BSS方法的目的是從觀測(cè)信號(hào)X(t)∈N×T中恢復(fù)源信號(hào)S(t)∈m×T，混合過程A∈N×m和源信號(hào)都是未知的[25]。其中：N是傳感器個(gè)數(shù)，T是采樣點(diǎn)數(shù)，m是源信號(hào)個(gè)數(shù)。不考慮測(cè)量噪聲，可以表示為

X(t)=AS(t)。

(1)

為了從觀測(cè)信號(hào)中分離出獨(dú)立源信號(hào)，需要尋找分離矩陣B∈m×N，使得

Y(t)=BX(t)=BAS(t)。

(2)

式中Y(t)=[y1(t),…,ym(t)]T∈m×T是源信號(hào)S(t)∈m×T的估計(jì)。如圖1所示為盲源混合和分離的過程。

1.2 EASI算法理論描述

傳統(tǒng)BSS方法為離線和批處理的算法，且只適用于平穩(wěn)環(huán)境。因此，研究者們提出了自適應(yīng)BSS方法，其中，Cardoso等[17]在1996年提出EASI算法，該算法基于連續(xù)更新的想法，采用相對(duì)梯度，通過迭代得到分離矩陣B(t+1)。

(3)

式中：I∈m×m為單位矩陣;λ(t)為第t時(shí)刻的步長;B(t+1)∈m×N取經(jīng)過多次迭代收斂的t+1時(shí)刻分離矩陣。

(4)

式中：V為特征X(t)的特征向量;U為對(duì)角矩陣，其對(duì)角線由X(t)的特征值組成。

傳統(tǒng)的EASI算法式(3)中的λ(t)是固定值。近年來，學(xué)者們提出了一系列的變步長方法，本文選取指數(shù)衰減型[26]和時(shí)間遞減型[27]。步長隨著時(shí)間t增加而減小，分別如式(4)和式(5)所示，使得EASI算法能夠兼顧穩(wěn)定性和收斂性：

(5)

(6)

1.3 基于EASI算法的線性時(shí)不變結(jié)構(gòu)工作模態(tài)參數(shù)識(shí)別

線性時(shí)不變結(jié)構(gòu)的質(zhì)量M、阻尼C和剛度矩陣K不隨時(shí)間改變，其動(dòng)力學(xué)方程為：

(7)

X(t)=ΦQ(t)。

(8)

2 基于MWVEASI算法的線性慢時(shí)變結(jié)構(gòu)工作模態(tài)參數(shù)在線識(shí)別

2.1 線性慢時(shí)變結(jié)構(gòu)在模態(tài)坐標(biāo)下的表示

根據(jù)動(dòng)力學(xué)理論，物理坐標(biāo)下線性時(shí)變結(jié)構(gòu)的動(dòng)力學(xué)方程為：

=F(t),t∈[tbegin,tend]。

(9)

式中：M(t),C(t)和K(t)∈N×N分別表示時(shí)變的質(zhì)量、阻尼和剛度矩陣?；凇皶r(shí)間凍結(jié)”和短時(shí)時(shí)不變理論，線性時(shí)變結(jié)構(gòu)在時(shí)間[tbegin,tend]內(nèi)，質(zhì)量、剛度和阻尼可視為是“凍結(jié)”(時(shí)不變)的，因而，式(9)可表示為一系列時(shí)不變結(jié)構(gòu)S′(ε)構(gòu)成的集合S′：

(10)

(11)

(12)

(13)

在每個(gè)窗內(nèi)，該時(shí)段內(nèi)的模態(tài)振型和固有頻率可作為最中間時(shí)刻的瞬時(shí)模態(tài)參數(shù)。

2.2 MWVEASI算法

滑窗(moving window)技術(shù)利用短時(shí)時(shí)不變理論，使得在每個(gè)窗內(nèi)非平穩(wěn)振動(dòng)響應(yīng)信號(hào)可以看作是平穩(wěn)的,線性時(shí)變結(jié)構(gòu)可以看成是時(shí)不變的。隨著時(shí)間的推移，窗口向右滑動(dòng)，增加新數(shù)據(jù)同時(shí)刪除舊數(shù)據(jù)，以此類推，如圖3所示描述了該過程。

將滑動(dòng)窗和自適應(yīng)BSS的變步長EASI算法相結(jié)合，提出滑動(dòng)窗變步長EASI(MWVEASI)的時(shí)變BSS算法。在每個(gè)窗內(nèi)，應(yīng)用變步長EASI方法進(jìn)行自適應(yīng)盲源分離。

2.3 基于變步長MWVEASI算法的線性時(shí)變結(jié)構(gòu)工作模態(tài)參數(shù)在線識(shí)別方法

在圖4的流程中，首先將響應(yīng)數(shù)據(jù)分成一系列數(shù)據(jù)窗。然后，在每個(gè)窗內(nèi)，每個(gè)時(shí)刻采用不同步長的EASI方法進(jìn)行迭代。將窗內(nèi)最后時(shí)刻的混合矩陣和源信號(hào)估計(jì)的結(jié)果看成是窗內(nèi)中間時(shí)刻的結(jié)果。接著，對(duì)EASI識(shí)別的源信號(hào)使用快速傅里葉變換(Fast Fourier Transform, FFT)或單自由度技術(shù)(如自然激勵(lì)法(Natural Excitation Technique, NExT)等)得到模態(tài)固有頻率和阻尼比。

2.4 方法的適用范圍和理論分析比較

基于滑動(dòng)窗變步長EASI的線性慢時(shí)變結(jié)構(gòu)工作模態(tài)參數(shù)識(shí)別方法的適用范圍如下：

(1)系統(tǒng)是線性慢時(shí)變的弱阻尼結(jié)構(gòu)；

(2)振動(dòng)響應(yīng)傳感器個(gè)數(shù)大于或等于可識(shí)別工作模態(tài)參數(shù)階數(shù)；

(3)振動(dòng)響應(yīng)信號(hào)的采樣頻率大于或等于可識(shí)別模態(tài)固有頻率的2倍；

(4)系統(tǒng)未知環(huán)境激勵(lì)應(yīng)該近似為平穩(wěn)白噪聲。

EASI方法只能識(shí)別線性時(shí)不變結(jié)構(gòu)的模態(tài)參數(shù)[19]。變步長的EASI方法通過可變的步長使得算法可以迅速收斂[20]。當(dāng)系統(tǒng)突變時(shí)，變步長EASI方法依舊可以收斂到平衡位置[22-23]。然而，一方面變步長EASI方法重新收斂需要時(shí)間，另一方面當(dāng)時(shí)變結(jié)構(gòu)參數(shù)隨時(shí)間持續(xù)變化時(shí)，EASI方法此時(shí)的步長不僅與時(shí)間相關(guān)，還應(yīng)與時(shí)變特性有關(guān)。只隨時(shí)間變化的步長將導(dǎo)致變步長EASI跟蹤時(shí)變的能力弱、識(shí)別精度低?；瑒?dòng)窗方法將當(dāng)前時(shí)刻數(shù)據(jù)和前后一定時(shí)延時(shí)刻的數(shù)據(jù)組成數(shù)據(jù)窗，在窗內(nèi)，基于“時(shí)間凍結(jié)理論”，數(shù)據(jù)非平穩(wěn)程度低，使得窗內(nèi)變步長EASI方法有較高的精度。滑動(dòng)窗過程加入了新數(shù)據(jù)并剔除舊數(shù)據(jù)，窗內(nèi)數(shù)據(jù)包含了當(dāng)前的時(shí)變信息，因而滑動(dòng)窗能夠跟蹤時(shí)變結(jié)構(gòu)的變化[24]。滑動(dòng)窗并不改變EASI算法本身的收斂性，因此滑窗變步長EASI收斂速度比滑動(dòng)窗固定步長EASI快。綜上所述，表1描述了EASI方法和滑動(dòng)窗方法相結(jié)合后的特點(diǎn)。

表1 基于滑動(dòng)窗變步長EASI與其他線性慢時(shí)變結(jié)構(gòu)工作模態(tài)參數(shù)識(shí)別的比較

3 仿真驗(yàn)證

3.1 仿真系統(tǒng)和數(shù)據(jù)集介紹

為驗(yàn)證所提出的MWVEASI算法的正確性和精度，在MATLAB/Simulink中，設(shè)計(jì)了一個(gè)線性時(shí)變的三自由度彈簧振子，如圖5所示。

圖5所示的三自由度彈簧結(jié)構(gòu)在慢時(shí)變條件下的動(dòng)力學(xué)方程為：

(14)

設(shè)置3個(gè)物塊的初始位移為零，阻尼比都為c1=c2=c3=0.01 N.s/m，剛度都是k1=k2=k3=1 000 N/m。質(zhì)量m2=m3=1 kg，m1隨時(shí)間緩慢變化如下：

(15)

m1(t)所受到的激勵(lì)F1(t)為高斯白噪聲，采樣間隔為0.025 s，采樣頻率為40 Hz，仿真時(shí)長為2 000 s。取50 s后的時(shí)變數(shù)據(jù)進(jìn)行研究，三自由度結(jié)構(gòu)在MATLAB/Simulink的建模如圖6所示。在Simulink模塊中利用Runge-Kutta算法求解得到位移響應(yīng)信號(hào)X(t)數(shù)據(jù)集，如圖7所示。更具體的仿真細(xì)節(jié)可見文獻(xiàn)[24]。最終得到的位移響應(yīng)信號(hào)數(shù)據(jù)集是一個(gè)3行80 001列的矩陣，即N=3，T=80 001。

3.2 評(píng)價(jià)指標(biāo)

理論的模態(tài)固有頻率、振型、阻尼比通過有限元仿真得到。采用模態(tài)置信度準(zhǔn)則(MAC)作為評(píng)價(jià)指標(biāo)衡量模態(tài)振型識(shí)別的精度，定義如下：

(16)

其中：φi為第i階理論值，φj為第j階識(shí)別值。0≤MAC≤1，MAC越接近1說明精度越高。

3.3 參數(shù)設(shè)置

滑動(dòng)窗的窗長L取1 024。變步長EASI方法的參數(shù)取值，即式(4)和式(5)分別為λ0=0.014,t0=200,α=0.0025,β=2。固定步長EASI方法取步長為0.001 4。B的初值設(shè)為：B0=0.6×rands(3,3)，其中rands()是MATLAB的一個(gè)函數(shù)，表示生成3行3列的每個(gè)元素是0～1之間隨機(jī)數(shù)的矩陣。ε=10-4，l_max=5000，t=l=i=1。

收斂和停機(jī)準(zhǔn)則：

(3)當(dāng)t=L時(shí)，i=i+1，進(jìn)入下一個(gè)窗，重置EASI方法初值，t=l=1，返回步驟(1)，直到i+L>T時(shí)停止。

3.4 工作模態(tài)參數(shù)識(shí)別結(jié)果

本文對(duì)兩種變步長EASI算法、相應(yīng)的兩種加滑動(dòng)窗變步長方法以及滑窗固定步長EASI的結(jié)果進(jìn)行了對(duì)比，圖8和圖9分別表示5種算法的頻率、MAC值識(shí)別結(jié)果與理論值的對(duì)比。表2統(tǒng)計(jì)了各種算法未識(shí)別的模態(tài)參數(shù)個(gè)數(shù)。

表2 未識(shí)別的模態(tài)參數(shù)百分比

3.5 識(shí)別結(jié)果分析

(1)從圖8b、圖8d、圖9b、圖9d可以看出，基于滑動(dòng)窗EASI方法能夠識(shí)別線性慢時(shí)變結(jié)構(gòu)的模態(tài)參數(shù)。而圖8a、圖8c、圖9a、圖9c則表明，不加窗的變步長EASI算法在識(shí)別時(shí)變快的第三階模態(tài)參數(shù)時(shí)效果不佳。對(duì)比圖8e與8a、圖8c，滑窗EASI方法時(shí)變快時(shí)識(shí)別的值仍有較高的精度，因而跟蹤性比變步長EASI好。

(2)對(duì)比圖8e與8b、圖8d，可以看出，滑窗固定步長EASI方法識(shí)別模態(tài)參數(shù)精度沒有滑窗變步長EASI方法高。從圖8b、圖8d、圖9b、圖9d可以看出，基于滑動(dòng)窗指數(shù)衰減變步長EASI算法識(shí)別精度比基于滑動(dòng)窗迭代遞減變步長EASI高。

(3)從圖8、圖9和表2可以看出，MWVEASI算法未能識(shí)別的模態(tài)參數(shù)時(shí)刻點(diǎn)的個(gè)數(shù)明顯少于MWEASI和VEASI，說明MWVEASI算法收斂速度快，更容易收斂、穩(wěn)定性更強(qiáng)。解決這些未識(shí)別的時(shí)刻點(diǎn)的辦法是將上一時(shí)刻點(diǎn)的分離矩陣收斂值作為該時(shí)刻點(diǎn)分離矩陣的迭代計(jì)算初值，以及改進(jìn)變步長格式和步長參數(shù)，加快算法收斂速度。

4 結(jié)束語

本文將滑動(dòng)窗方法與變步長EASI方法相結(jié)合，應(yīng)用于識(shí)別線性慢時(shí)變結(jié)構(gòu)的工作模態(tài)參數(shù)。在質(zhì)量慢時(shí)變的三自由度彈簧振子結(jié)構(gòu)的仿真數(shù)據(jù)集上的識(shí)別結(jié)果表明，兩種不同的MWVEASI方法都能夠很好地識(shí)別線性慢時(shí)變結(jié)構(gòu)的模態(tài)振型、固有頻率，與滑窗固定步長EASI、變步長EASI算法相比，MWVEASI能夠更好地跟蹤時(shí)變特性、識(shí)別精度高、收斂速度快。

但是，變步長也存在方式和參數(shù)選擇的問題；滑動(dòng)窗方法也帶來了計(jì)算量的增加和窗長設(shè)置的問題；如何實(shí)現(xiàn)窗與窗之間的遞推，減少時(shí)空復(fù)雜度也是進(jìn)一步的研究方向。與此同時(shí)，通過具體的工程實(shí)例的實(shí)驗(yàn)數(shù)據(jù)來驗(yàn)證所提出的新方法也是進(jìn)一步的工作方向。