變運營環境下基于混合主成分分析的結構損傷識別方法

黃海賓，臧敬剛

(1. 河北工業大學 土木與交通學院，天津 300401； 2. 河北工業大學 河北省土木工程技術研究中心，天津 300401)

0 引 言

近年來，基于振動的損傷識別方法在結構健康監測領域得到了廣泛研究[1-2]。該類方法的基本依據是：損傷會引起結構動力特性(如頻率、振型等模態參數)的變化，故可將動力特性作為結構的損傷特征，并通過動力特性變化判斷其損傷狀態[3-5]。然而，服役期內的運營環境(如溫度、濕度、風速等因素)變化同樣會引起結構動力特性的變化，這可能會掩蓋損傷造成的動力特性變化，導致損傷識別效果不理想[6-11]。因此，在結構損傷識別中，剔除運營環境變化對損傷特征的影響至關重要。

從損傷特征中剔除運營環境變化的影響主要有兩類方法[12]：①顯式方法，建立環境變量與損傷特征之間的關系模型，進而量化運營環境變化的影響并予以剔除；②隱式方法，將環境變量視作隱藏變量嵌入到損傷特征中，進而估計運營環境變化的影響并予以剔除。隱式方法因無需測量環境變量卻能考慮運營環境變化的影響，進而凸顯出了實用價值，其中以主成分分析(principal-component analysis, PCA)應用最廣。A. M. YAN等[6]將固有頻率作為損傷特征，使用PCA成功剔除了環境效應對其損傷識別過程的影響；G. COMANDUCCI等[7]將PCA應用于某懸索橋，有效剔除了風速對固有頻率的影響；F. UBERTINI等[8]提出將PCA與多元線性回歸相結合，建立統計模型以剔除運營環境變化的影響，并采用某鐘樓的固有頻率數據進行了有效性驗證；A. I. OZDAGLI等[9]使用振型和固有頻率作為損傷特征，結合PCA剔除了溫度變化的影響，提升了對某三層框架的損傷識別能力。

盡管基于PCA的損傷識別方法應用眾多，但其適用范圍僅局限于數據近似滿足高斯分布且線性相關。目前，針對PCA在非線性問題中的局限性，學界已提出相關改進方法，A. M. YAN等[10]將PCA擴展為局部PCA，對非線性相關性進行分段線性化，剔除了運營環境變化對固有頻率的非線性影響；E. REYNDERS等[11]針對變運營環境的非線性影響，提出了改進的核PCA方法，并采用某預應力混凝土橋的固有頻率數據對其效果進行了驗證。此外，非高斯分布在損傷特征數據中也普遍存在，趙人達等[13]和薛剛等[14]分別研究了不同太陽輻射強度下的混凝土箱梁溫度場分布特性，結果表明箱梁底板、腹板和頂板的溫度均呈非高斯分布，這會進一步造成損傷特征數據的非高斯分布。在結構損傷識別中，綜合考慮數據的非高斯分布和非線性相關十分關鍵。

為有效處理運營環境變化下損傷特征數據的非高斯分布和非線性相關等問題，筆者提出一種基于混合PCA的結構損傷識別方法。利用高斯混合模型(gaussian mixture model, GMM)對存在非高斯分布和非線性相關的多維損傷特征數據進行建模，將其聯合概率密度函數表示為多個局部高斯分量的線性組合，使得每個局部高斯分量內的數據之間滿足線性相關性；然后，對所有局部高斯分量分別建立相應的PCA模型；再計算所有PCA模型殘差部分的馬氏平方距離和歐氏平方距離，將樣本屬于各高斯分量的后驗概率作為權系數，對兩類距離分別進行加權標準化以求得結構的綜合損傷指標。最終，通過質量彈簧系統仿真和木桁架橋試驗對其有效性進行驗證。

1 傳統PCA損傷識別

作為一種常用的多維數據分析方法，PCA可對由多個互相關特征變量構成的原始數據集進行線性變換以提取主成分，在保留絕大部分原始信息的前提下，降低變換后數據集的特征維度。原始數據集中各特征變量間的相關性越強，PCA的降維效果越顯著。

令X=[x1,x2,…,xN]表示損傷特征矩陣，任意樣本向量x均包含m個特征變量，共有N個樣本。采用結構的多階頻率作為損傷特征，且各階頻率間存在相關性，故可通過PCA對特征數據集進行建模。PCA建模是利用如下特征值分解實現[15-16]：

∑=E[(x-μ)(x-μ)T]=QΛQT

(1)

式中：μ為數據集的均值向量；∑為數據集的協方差矩陣；Λ=diag(λ1,λ2,…,λm)為所有降序排列特征值構成的對角矩陣；Q=[q1,q2,…,qm]為與特征值對應的特征向量構成的標準正交矩陣；E(·)為期望算子。

PCA模型中特征值λi越大，則表明第i個主成分所包含的原始數據信息越多。當特征變量間存在強相關性時，僅需前面少數幾個主成分即可近似重構原始數據。通常情況下，主成分個數d的選取原則如下：

(2)

式中：r為百分率，一般取95%，即主成分貢獻率宜至少達到95%。當已知影響損傷特征的主導環境因素的個數時，可將d取為主導環境因素的個數。

(3)

(4)

(5)

(6)

當數據服從或近似服從高斯分布時，令α表示顯著性水平，指標DMah的閾值TMah由下式確定[17]：

(7)

(8)

2 混合PCA損傷識別

2.1 損傷特征的聯合概率密度擬合

GMM利用多個高斯分量的線性組合對多維數據的概率密度函數進行擬合[18]。理論上講，如果高斯分量的數目足夠多，GMM能夠準確擬合任意分布數據的概率密度函數。

經GMM擬合，多維損傷特征數據的概率密度函數具有如下形式：

(9)

(10)

(11)

進一步，可將GMM的最優參數估計問題轉化為如下優化問題：

(12)

采用期望最大(expectation-maximization, EM)算法[19]求解該問題，首先設定GMM的初始參數為Θ(0)，然后分E步驟和M步驟進行迭代計算。E步驟依據當前參數Θ(i)，計算第j個樣本向量xj屬于第k個高斯分量的后驗概率：

(13)

M步驟迭代更新參數Θ(i+1)：

(14)

(15)

(16)

重復E步驟和M步驟，直至收斂，即可得到GMM的最優參數。

在GMM擬合中，確定最佳的高斯分量數目是關鍵所在，采用可綜合權衡模型擬合度與復雜度的貝葉斯信息準則對高斯分量的數目進行選擇。

2.2 混合PCA建模與損傷指標計算

多維損傷特征數據經GMM擬合后，可獲得所有高斯分量的均值向量和協方差矩陣，且任一高斯分量中的損傷特征之間存在線性相關性。由式(1)可知：PCA建模的實質是對協方差矩陣的特征值分解，故可依次對第k個高斯分量的協方差矩陣∑k分別進行特征值分解，從而得到K個PCA模型及其對應的用于計算損傷指標的殘差子空間。

由于每個高斯分量對應的殘差子空間不同，以其計算而得的損傷指標及閾值也會存在差異，故分別對馬氏平方距離和歐氏平方距離進行加權標準化，作為混合PCA框架下的損傷指標。對于任意樣本向量x而言，首先計算第k個高斯分量對應的損傷指標，即馬氏平方距離DMah,k和歐氏平方距離DEuc,k；其次，通過相應的閾值TMah,k和TEuc,k分別對損傷指標進行標準化處理，即將損傷指標除以其閾值；最后，對所有標準化后的損傷指標進行加權求和，其權系數為該樣本來自GMM各高斯分量的后驗概率p(k|x)。具體表達式為：

(17)

(18)

由于進行了加權標準化處理，兩個綜合損傷指標的閾值均為1。當損傷指標大于1時，判斷結構處于損傷狀態；反之，則判斷結構處于無損狀態。

2.3 損傷識別步驟

基于混合PCA的損傷識別方法，其實施過程分為離線建模和在線監測兩個階段。

2.3.1 離線建模

離線建模階段，利用無損結構的損傷特征數據集進行混合PCA建模，步驟如下：

1)采用貝葉斯信息準則確定GMM的最佳高斯分量數目；

2)通過GMM對損傷特征數據的概率密度函數進行擬合；

3)分別對GMM中各高斯分量建立PCA模型并獲得其對應的殘差子空間；

4)計算GMM中各高斯分量所對應的馬氏平方距離和歐氏平方距離的閾值。

2.3.2 在線監測

在線監測階段，計算混合PCA框架下當前樣本的損傷指標以判斷結構損傷狀態，步驟如下：

1)計算當前樣本對應于GMM中各高斯分量的馬氏平方距離和歐氏平方距離；

2)計算當前樣本來自于GMM中各高斯分量的后驗概率；

3)計算當前樣本的加權標準馬氏平方距離和加權標準歐氏平方距離作為損傷指標；

4)判斷損傷狀態，當損傷指標大于1時，判斷結構處于損傷狀態，反之則判斷結構處于無損狀態。

3 數值仿真

利用四自由度質量彈簧系統生成仿真數據，用以驗證所提方法在變運營環境下的損傷識別能力。

質量彈簧系統共由4個質量塊和5根彈簧組成，如圖1。質量塊的質量為：m1=m2=m3=m4=2 kg；彈簧的剛度與溫度呈分段線性關系：

圖1 質量彈簧系統

(19)

(20)

式中：ki,i=1,2,…,5為彈簧剛度；T為溫度。

在0 ℃兩側，所有彈簧的剛度隨溫度的線性變化規律不同；第3根彈簧的剛度隨溫度的線性變化規律與其它彈簧的變化規律均不一致，這可用于模擬非線性效應。

以某橋梁監測系統采集的為期1年的空氣溫度數據(采樣頻率1 Hz，1 h內會得到多個樣本，對這些樣本取平均值，得到8 760個平均樣本)為基準，等比例伸縮至[-20 ℃, 40 ℃]區間，作為數值仿真的溫度輸入值。通過無放回隨機取樣，首先取出5 400個溫度樣本輸入質量彈簧系統，計算全部4階固有頻率即fn1、fn2、fn3和fn4，疊加一定程度的高斯噪聲后，作為訓練數據集；然后，分7次依次各取出480個溫度樣本輸入質量彈簧系統，且每次均對第2根彈簧的剛度進行一定程度折減，在模擬運營環境變化和損傷同時存在情況下，計算全部4階固有頻率并疊加與訓練數據集程度相同的高斯噪聲后，作為測試數據集(共計7種工況)，如表1。

表1 測試數據集工況

該系統損傷前后的固有頻率變化如圖2。經對比分析可知：運營環境變化引起的固有頻率變化會大部分甚至完全掩蓋損傷引起的固有頻率變化，故僅通過頻率變化難以識別運營環境變化下的結構損傷。

圖2 質量彈簧系統的固有頻率變化

由貝葉斯信息準則確定最佳高斯分量數目為4，利用GMM對訓練數據集的概率密度函數進行擬合(以fn1和fn2為例)。由圖2可知：無損狀態下fn1和fn2的變化幅度較大。圖3(a)和圖3(b)分別為fn1和fn2的頻率直方圖與GMM擬合概率密度函數間的對比，可知GMM能有效擬合非高斯分布數據的概率密度函數；圖3(c)為fn1和fn2間的散點圖，其表現出非線性相關性；圖3(d)為GMM擬合的二維概率密度函數等高線圖，可知GMM亦能有效處理數據中的非線性問題。

圖3 頻率fn1、fn2及 fn1與fn2間的GMM擬合效果

圖4 混合PCA的損傷識別結果

表2 PCA與混合PCA的損傷識別率

由表2可知：對工況Cn0，混合PCA的兩種損傷指標對應的虛警率(對正常結構識別出損傷)較傳統PCA均略低；對工況Cn1，混合PCA的識別效果優于傳統PCA，但此時的微小損傷引起的固有頻率變化不明顯，漏警率(對損傷結構未識別出損傷)較高；對工況Cn2，混合PCA的損傷識別能力較傳統PCA有了顯著提升；對工況Cn3～Cn6，混合PCA的損傷識別率均接近或超過90%，且均大于傳統PCA的損傷識別率，足見筆者方法的優越性。

4 試驗驗證

采用由芬蘭的Kullaa課題組所完成的木桁架橋試驗[20]進一步驗證混合PCA的損傷識別效果。

橋總質量為36 kg，采用橡膠支座支撐，如圖5。對結構施加隨機白噪聲激勵，并在不同位置安裝15個加速度計，以測量結構的加速度響應，采樣頻率為256 Hz，每次測量時長為32 s。測量期間，實驗室內的溫度和濕度不斷發生變化。

圖5 木桁架橋的試驗布置

利用隨機子空間方法對每次測量的振動數據進行模態參數識別，并采用較為可靠的第6、7、8、10、12和13階固有頻率(為后續方便，將它們依次記為ft1、ft2、ft3、ft4、ft5和ft6)作為損傷特征。在無損狀態下，共獲得1 871個有效的固有頻率樣本，通過無放回隨機取樣，選取1 671個樣本作為訓練數據集，剩余200個樣本作為測試數據集的一部分；在距中跨左側600 mm處的結構頂部附加不同質量，以模擬不同的損傷程度，共獲得105個有效的固有頻率樣本，作為測試數據集的另一部分。因此，測試數據集中共包含6種工況，如表3。

表3 測試數據集工況

木桁架橋在環境影響下的固有頻率變化如圖6，圖中虛線前為無損狀態，虛線后為損傷狀態。經對比計算可知：由環境變化引起的固有頻率變化完全掩蓋了損傷引起的固有頻率變化，故僅通過頻率變化無法識別運營環境變化下的結構損傷。

由貝葉斯信息準則確定最佳高斯分量數目為8，利用GMM對訓練數據集的概率密度函數進行擬合。以ft2和ft6為例，由圖6可知兩頻率在無損狀態下的變化幅度較大；圖7(a)和圖7(b)分別為ft2和ft6的頻率直方圖與GMM擬合概率密度函數間的對比，可知GMM能有效擬合非高斯分布數據的概率密度函數；圖7(c)為ft2和ft6間的散點圖，其表現出非線性相關性；圖7(d)為GMM擬合的二維概率密度函數等高線圖，可知GMM能有效處理數據中的非線性問題。

圖6 木桁架橋的固有頻率變化

圖7 頻率ft2、ft6及ft2與ft6間的GMM擬合效果

圖8 混合PCA的損傷識別結果

表4 PCA與混合PCA的損傷識別率

由表4可知：傳統PCA中的歐氏平方距離并沒有識別出結構損傷；傳統PCA中的馬氏平方距離雖有一定的識別效果，但并沒有識別出工況Ct3的損傷；混合PCA中的加權標準歐氏平方距離的識別效果較傳統PCA有了一定改善，但其對工況Ct3的識別效果仍不理想；混合PCA中的加權標準馬氏平方距離對工況Ct2識別率達到了60.87%，對工況Ct3～Ct5識別率則均達到了100%，識別結果最好。綜上，由于可有效處理非高斯分布和非線性相關等問題，混合PCA在損傷識別能力方面較傳統PCA有顯著提升。

7 結 論

筆者提出了變運營環境下基于混合PCA的結構損傷識別方法，并通過質量彈簧系統仿真和木桁架橋試驗進行了有效性驗證。得出如下結論：

1) 運營環境變化會引起結構損傷特征的非高斯分布和非線性相關等問題，致使傳統PCA的損傷識別效果較差。

2) 在合理的高斯分量數目下，GMM能準確擬合多維損傷特征數據的聯合概率密度函數，解決非高斯分布和非線性相關等問題。

3) 混合PCA的損傷識別能力較傳統PCA而言有顯著提升，其中以加權標準馬氏平方距離的損傷識別能力最為優異。