基于密度峰值聚類算法的模態參數識別

王飛宇，胡志祥，黃 瀟

(1.合肥工業大學 土木與水利工程學院，合肥 230009； 2.中國電子科技集團公司第三十八研究所，合肥 230088)

僅通過系統的響應信號識別模態參數一直是個有挑戰性的研究。盲源分離法(Blind Source Separation, BSS)已經被廣泛的應用于模態分析領域。Kerschen等首次將單模態信號作為BSS算法中的源信號，并利用獨立分量分析(Independent Component Analysis，ICA)提取多個單模態信號，得到結構每階模態振型。與傳統時域方法相比，BSS具有無需先驗信息、計算簡單、非參數化等優點。隨后，一些學者提出了新的盲源分離方法都得到較好的結果，例如：二階盲辨識方法(Second Order Blind Identification, SOBI)、多源提取方法(Algorithm for Multiple Unkown Signals Extraction, AMUSE)等。

然而，文獻[1-7]中的算法都只能解決BSS算法中的正定問題，即需要滿足傳感器數目(m)大于或等于結構自由度數目(n)。當m

本文首先利用短時傅里葉變換將信號從時域轉換到時頻域。針對源信號數量未知和結構高階模態混疊的問題，利用2014年發表于science的聚類算法識別模態振型[13]，它是一種基于密度峰值的聚類方法(Desity Peaks Clusering Algorithm, DPCA)。DPCA算法可以在源信號數量未知的情況下，通過決策圖直觀得到聚類數目，并且有效忽略頻率混疊點和噪聲點對聚類精度的影響。最后，利用SL0算法重構每階模態的時頻域信號，通過逆短時傅里葉變換得到單模態時域信號，進而計算出模態頻率。

1 盲源分離與模態參數識別

考查線性時不變結構系統，其在激振力作用下的微分方程為

(1)

式中：M，C，K分別為n×n階質量，阻尼和剛度矩陣；f(t)為激振力；x(t)為位移響應。將x(t)轉化到模態坐標形式為

x(t)=Φq(t)

(2)

式中：Φ為模態振型矩陣，其每一列向量代表結構的某一階振型；q(t)為模態坐標振動信號。各階模態坐標振動信號可表示為

qi(t)=Bie-ξiωn,itsin(ωd,it+φi)

(3)

式中：ωn,i和ωd,i分別為無阻尼和有阻尼的情況下每階振型的固有頻率；B和φ為兩個由初始條件確定的常數；ζi為第i階振型的阻尼比。僅通過系統的響應信號識別模態參數，即僅利用振動響應x(t)提取結構模態振型、頻率和阻尼比。

上述模態參數識別問題可以與盲源分離問題對應起來。盲源分離是指在源信號與混合通道參數均未知的情況下，僅利用傳感器測量信號計算各源信號的信號處理方法。BSS的數學模型為

X(t)=AS(t)

(4)

式中：X(t)∈Rm為觀測信號，m為觀測信號數量；A∈Rm×n為未知列滿秩矩陣；S(t)為未知n維源信號。對于正定問題，即m=n，盲源分離的主要任務是僅通過觀測信號找到混合矩陣A并求A-1，然后由式(1)解出S(t)。

對比式(2)和式(4)可以看出，盲源分離的數學模型與模態分析時域模型有相似之處。由于源信號和混合矩陣的先驗信息未知，兩者都是僅通過結構系統的輸出信號求解方程。模態響應q(t)相當于式(4)中的S(t)，而位移響應x(t)相當于式(4)中的X(t)。本文所要解決的為欠定盲源分離問題，即m

2 模態振型矩陣計算

首先利用調頻信號算例說明模態振型矩陣計算原理。設混合信號由兩個簡諧信號和一個調頻信號組成，簡諧信號q1和q2頻率分別為10 Hz和20 Hz，調頻信號q3從5 Hz增加到15 Hz，調頻信號的頻率范圍保證了信號在時頻域有混疊成分；采樣頻率100 Hz；加入混合信號能量10%的白噪聲。設矩陣Φ為

(5)

利用短時傅里葉變換將混合信號x(t)轉換到時頻域上，得到稀疏信號

x(t,f)=Φq(t,f)

(6)

該混合信號的時頻圖如圖1所示。式(6)中源信號數量為3，觀測信號數量為2。式(6)亦可寫作

(7)

圖1 混合信號時頻圖Fig.1 The mixing signals in the time-frequency domain

假設存在一點(t1,f1)，僅存在源信號q1，那么式(7)可以簡化為式(8)和式(9)

Re(x1(t1,f1))=a11Re(q1(t1,f1))

(8)

Re(x2(t1,f1))=a21Re(q1(t1,f1))

(9)

式中：Re為取實部，由式(8)和式(9)可知

Re(x1(t1,f1))/Re(x2(t1,f1))=a11/a21

(10)

由式(10)可知，通過混合信號時頻點實幅值的比值可得到混合矩陣第一列向量，重復這一步驟可以得到混合矩陣每一列向量。式(10)中選取虛部和實部的結果是一樣的[14]。繪出兩個混合信號在時頻域上的實部幅值散點圖，如圖2所示。圖中散點的排列方向趨近于Φ矩陣各列向量的方向，如圖虛線所示。通過聚類算法求出每一組的聚類中心，即為振型矩陣的每一列向量。

圖2 混合信號幅值散點圖Fig.2 Scatter plot of the amplitude of the mixing signals

圖2中，位于虛線上的點稱為單源信號點，而遠離虛線的點稱為多源信號點。多源信號點會影響聚類算法的精度。為提高聚類算法精度與效率，宜先剔除多源信號點。當同一時頻點對應的實部向量與虛部向量的夾角低于某一閾值時，可以認為該點為單源信號點[15]。實部，虛部和夾角滿足下列關系式

(11)

式中：Im為虛部；θ為混合信號實部向量與虛部向量的夾角閾值，2°即可滿足要求。

選取滿足式(11)的時頻點可以保證每個源信號互相獨立，在時頻域上無混疊成分。圖3為去除多源信號點的混合信號幅值散點圖。圖中三類標記的數據點分為三個排列方向，即混合矩陣的列向量。可以看出，通過設置實部與虛部角度的閾值可以很大程度上減少計算量并且提高聚類精度。對這些時頻點進行聚類可計算出模態振型矩陣，聚類方法在下一節進行介紹。

圖3 去除多源點混合信號幅值散點圖Fig.3 Scatter plot of the amplitude of the mixing signals after eliminating source points

3 DPCA算法

本文利用一種基于峰值密度的聚類算法估計模態振型，該算法應用于模態分析的方面有兩個優點：①不需預先知道結構自由度數或活躍模態數，通過決策圖直觀選出聚類中心和聚類數目；②可以有效地分離噪聲點，從而保證模態振型計算的精度。

算法的核心思想是尋找被低密度區域分隔的高密度峰值點。對于一個數據集，聚類中心總是被低局部密度的數據點包圍。DPCA主要分為2步：①通過計算每個數據點的局部密度ρi和與更高密度點之間的距離δi，確定聚類中心；②將所有數據點劃分到對應的類簇中。

3.1 確定聚類中心

局部密度ρi的定義為

ρi=∑jχ(dij-dc)

(12)

其中，

(13)

式中：dc稱為截斷距離，為自定義參數。所選dc使得每個數據點的平均局部密度為數據點總數的1%～2%。本文所用距離公式為1-|cosα|，α為數據點對應的向量夾角。

對于數據點i，若存在更高密度的數據點，δi由式(14)計算

δi=min(dij),j∈Is

(14)

式中：Is為所有局部密度大于數據點i的序號的集合。即，δi表示在所有局部密度大于xi的數據點中，與xi距離最小的數據點與xi之間的距離。當數據點i有最大局部密度時，Is=?，δi取最遠數據點與xi之間的距離。

因此，每一個數據點xi有對應局部密度和距離值：(ρi,δi)。在二維圖中畫出每一組(ρi,δi)，即為決策圖。當xi同時有較大的ρ值和δ值時，則被認為是數據集的聚類中心。考慮圖3中的算例，將所有數據點對應的(ρi,δi)在平面上表示出來(ρ為橫軸，δ為縱軸，ρ和δ已被歸一化)，如圖4所示。

圖4 圖3中數據點決策圖Fig.4 The decision graph for the data inFigure 3

圖4中有3個標記點“脫穎而出”，同時擁有較大的ρ值和δ值，表示源信號有3個。每一種標記類型點對應著圖2中的一個排列方向。由圖4可知，正方形聚類中心有最大局部密度ρ，此時式(14)中的Is為空集。因此，正方形聚類中心的δ取最遠的數據點到它的距離。三角形與星形聚類中心不具有最大的局部密度，但局部密度比它們大的點屬于其他類簇，因為聚類中心在該類簇中有最大的ρ值。因此，3種標記點都有較高的δ值。對于圓點，比它ρ更大且距離最近的點與其屬于同一類簇，所以圓點δ值較小。

若決策圖中無法直觀地得到聚類數目時，可通過計算γ值確定類簇數目

γi=ρiδi,i∈Is

(15)

圖5 數據點γ值分布圖Fig.5 The value of γ in decreasing order

3.2 剩余數據點歸類

考慮2.1中的算例，在確定了3個聚類中心后，再對圖3中所有的圓點數據進行歸類。決策圖(見圖4)中各聚類中心ρ值的大小關系為：三角形的密度<星形的密度<正方形密度。對比圖3與圖4也可發現，在正方形類簇對應的方向上，數據點聚集得最為密集，而三角形與星形點略為稀少。

4 SL0算法基本原理

得到局部振型矩陣后，式(6)變為一個欠定的線性方程組。由于源信號在時頻域具有稀疏性，附加一定限制條件可以尋找最稀疏解。于是求解式(6)轉化為一個最優化問題，對于時頻點(t,f)，求解模型為

(16)

l0范數最小優化方法是求解該問題的方法之一。但是由于l0范數對噪聲極為敏感，任何大于零的數都會被當作非零數。通常用l1范數代替l0范數，轉化為一個凸優化問題，可以用線性規劃方法進行求解，但該方法復雜度高，計算量大。因此Mohimani等通過最速下降法和梯度投影原理提出用平滑函數趨近信號的l0范數，把最優化問題轉化為平滑函數的極值問題，稱為平滑l0范數重構算法。它可以在保證同樣精度的情況下，重構速度比基追蹤算法[16]快2～3倍。

SL0的核心思想是用一個連續光滑函數去估計不連續的l0范數。此光滑函數采用均值為0的高斯函數

(17)

式中：s為復數；σ為函數的光滑程度，數值越大越函數光滑。當σ趨近零時

(18)

l0范數計算向量中的非零元素，而SL0是計算向量中零元素的數量。對于有n維的任意向量s，其近似連續函數可以寫成

(19)

最后l0范數最小值問題歸為連續函數Fσ(s)最大值問題。

‖s‖0≈n-Fσ(sj)

(20)

對于上文所述算例，在時頻點(t,f)處，觀測量有兩個，未知量有3個。SL0算法即是要找出一組q1(t,f),q2(t,f),q3(t,f)，在滿足式(6)的同時使其平滑l0范數最小(或‖s‖0最大)。

5 模態識別算例

采用一個6層剪切框架模型作為算例，各層剛度和質量分別為176.729 MN/m和100 t；第一階阻尼比為0.05；采樣頻率取50 Hz，采樣時間取40 s；各層受白噪聲激勵，采用龍格-庫塔法進行振動響應計算；位移傳感器數量取4個，少于結構自由度數量，分別布置在2層、4層、5層、6層。對響應信號進行短時傅里葉變換。為了減小聚類計算量，先利用式(11)去除在時頻域發生頻率混疊的時頻點，再進行后續識別。

5.1 模態振型計算

在計算結構陣型時，選取dc使得每個數據點的平均鄰居個數為數據點總數的1%。圖6和圖7分別為峰值聚類算法的決策圖與γ值降序分布圖。

圖6 仿真算例決策圖Fig.6 The decision graph of the simulation example

圖7 仿真算例γ值分布圖Fig.7 The value of γ in decreasing order for the simulation example

圖6中有6個三角形點同時擁有較大的ρ值和δ值，表示數據點被分為6類簇。圖7也可以看出，有6個點γ值較大。可判定聚類中心為6個。

由于環境激勵和阻尼會使散點圖中有很多噪聲點，使高階模態難以識別。為了區分噪聲點，DPCA方法中引入了邊界區域的概念。邊界區域的是屬于某一類簇但與屬于其它類簇的數據點距離不超過dc的數據點的集合。在邊界區域中，定義最大的局部密度點對應的局部密度為ρb。若類簇中數據點的密度值高于ρb，則視為核心點；剩下的則被視為類簇光暈點。將圖6中三角形點作為聚類中心，得到的核心點與光暈點如表1所示。例如表1中的類簇2，屬于該類簇的點數共有1 002個，其中大于邊界密度的核心點有255個，小于邊界密度的光暈點共有747個。取核心點平均方向作為類簇方向，可有效降低模態振型識別誤差。

為表明DPCA聚類方法的優點，用模態置信準則(Modal Assurance Criterion, MAC)評估模態振型識別結果，對比了三種方法聚類識別出的振型結果，各階振型MAC如表2所示。

由表2可知，K均值和分層聚類識別的第一階和第二階振型并不精確。K均值聚類初始值選取不佳可能導致算法無法收斂到最優解，且需要預先確定類簇數目；分層聚類雖然不需要輸入類簇數目，但是一旦聚類被分裂或合并后，就無法返回，聚類質量受限制[17]。而DPCA聚類可以把數據點分為核心點和光暈點，光暈點中包含噪聲點。在計算聚類均值時，僅采用核心點，保證了計算模態振型的精度。識別振型與理論振型的對比如圖8所示，可見各階模態振型均能被精確識別。

表1 DPCA聚類后數據點分布情況Tab.1 The distribution of the points

表2 三種聚類方法計算振型MAC值Tab.2 The MAC obtained by three clustering algorithms

圖8 模態振型識別值與理論值對比Fig.8 The mode shapes of the identification value and the theoretical value

為了驗證明該算法的抗噪性，對傳感器信號加入白噪聲。信噪比分別取20 dB，10 dB，5 dB，2 dB，不同信噪比下的DPCA決策圖如圖9所示。

從圖9可以看出，隨著傳感器信號信噪比減小，決策圖中處于“特殊位置”的標記點數目沒有發生變化，始終是6個。將圖中標記點作為聚類中心所識別的模態振型MAC值，如圖10所示。

圖10中各階振型的MAC值均在0.99～1，說明該算法識別模態振型結果受噪聲影響較小。且圖9和圖10證明了PDCA在較高噪聲水平下仍可以準確選擇出結構的自由度數目和識別模態振型。

圖9 各種信噪比下的決策圖Fig.9 The decision graph under various SNR

圖10 各種信噪比下識別的模態振型Fig.10 The mode shapes MAC under various SNR

5.2 單模態信號重構

得到模態振型后，采用SL0算法重構各階單模態信號。重構后的各階單模態振動信號時頻域能量分布如圖11所示。

圖11中最大能量對應的頻率值由直線標出。為提高頻域分辨率，可通過逆短時傅里葉變換將時頻域信號轉換到時域，再由快速傅里葉變換得到更精確的結果。各階模態坐標信號如圖12所示，對應傅里葉變換如圖13所示。識別出的各階模態頻率見表3。由識別結果可見各階模態頻率識別精度較高。

表4給出了不同信噪比下算法識別出的模態頻率包含的誤差，各信噪比下模態頻率識別結果較穩定，誤差基本在1%左右，只有個別誤差超過2%，這說明算法對噪聲有較好的魯棒性。

除上述數值算例外，利用本文算法，也精確地識別出了ASCE benchmark模型[18]的模態振型和頻率，限于文章篇幅較長，未在文章中展示。

表3 頻率識別值與理論值Tab.3 The mode frequency of the identification value and the theoretical value

表4 不同信噪比時的識別模態頻率識別誤差Tab.4 Modal frequency error under various SNRs

圖11 各階模態信號時頻域三維圖Fig.11 Time-frequency representation of recovered modal coordinate signals

圖12 各階模態時域信號Fig.12 Time domain represent of the recovered modal responses

圖13 各階模態頻域信號Fig.13 Time frequency represent of the recovered modal responses

6 結 論

本文利用基于峰值密度的聚類算法識別模態振型，然后通過可以快速重構稀疏信號的SL0算法重構各階單模態振動信號，并識別各階模態頻率。在數值算例中，DPCA算法可以從決策圖和γ值分布圖中直觀地識別出源信號數量。在計算聚類中心時，DPCA算法將所有數據點分為核心點與光暈點，然后從核心點中得到更精確的聚類均值。

在重構信號時，充分利用了結構振動響應信號在時頻域上的稀疏特性，采用SL0算法，精確高效地重構出單模態振動信號。數值算例顯示，在較低信噪比下，算法仍可以準確識別結構模態振型和模態頻率，對噪聲有較好的魯棒性。算法適用于傳感器數量小于結構自由度數的盲源分離問題，為結構模態參數識別提供了新的選擇。