鄭 泓，段忠東

（哈爾濱工業(yè)大學（深圳）土木與環(huán)境工程學院，廣東深圳518055）

引言

土木工程結構在服役期間容易受到環(huán)境腐蝕、車輛荷載、自然災害等的共同作用而出現(xiàn)功能退化，嚴重時甚至會導致結構失效而引發(fā)災難性事故。為保證人員生命和財產(chǎn)安全，有必要對在役結構的狀態(tài)進行監(jiān)測［1-2］。近幾十年來，基于振動的結構損傷預警方法已成為結構健康監(jiān)測（Structural Health Monitoring，SHM）的重要手段［3-5］。該方法認為損傷會改變結構動力參數(shù)（如模態(tài)頻率等），通過比較損傷前后這些參數(shù)的變化可以識別損傷。但是，環(huán)境因素（如溫度、車載等）的改變也會引起結構動力參數(shù)的變化［6-9］，這種變化往往會掩蓋結構的真實損傷。如果不消除損傷預警過程中環(huán)境因素的影響，SHM的可信度將大打折扣。為此，考慮環(huán)境影響的損傷預警方法開始受到學者的關注。

國內(nèi)外學者對這一問題開展了大量研究。Sohn和Farrar［10］將 消 除 環(huán) 境 因 素 影 響 的 過 程 稱 為“數(shù)據(jù)標準化（Data Normalization）”。根據(jù)是否需要測量環(huán)境信息，該過程分為基于模型的方法和非模型的方法。前者需要借助采集的環(huán)境信息建立結構動力參數(shù)與環(huán)境變量之間相互關聯(lián)的數(shù)學模型，以量化環(huán)境因素對結構動力參數(shù)的影響。Peeters和De Roeck［6］采 用ARX模 型 建 立 了Z24橋 的 環(huán) 境 溫度與模態(tài)頻率的關系；文獻［11-12］先后采用支持向量回歸模型（Support Vector Regression，SVR）和人工神經(jīng)網(wǎng)絡（Artificial Neural Network，ANN）量化溫度對香港汀九橋模態(tài)頻率的影響。Yang等［13］采用線性回歸方法建立潤揚長江大橋的模態(tài)頻率和溫度的回歸關系?；谀Ｐ偷姆椒槔斫猸h(huán)境因素如何影響結構物理特性的機理奠定了基礎，但實際工程中某些環(huán)境變量很難測量，準確建立環(huán)境因素和結構動力參數(shù)之間的關系難度很大。非模型的方法無需利用環(huán)境信息，而是將環(huán)境變量當作隱藏變量，通過機器學習的手段從大量樣本中找出環(huán)境因素影響的內(nèi)在規(guī)律［14］?；谥鞒煞址治觯≒rinciple Component Analysis，PCA）的方法［15］假定環(huán)境因素改變是結構動力參數(shù)出現(xiàn)變化的主要原因，采用PCA可以將環(huán)境因素以主成分的形式分離出來。林友新等［16］認為結構模態(tài)頻率受環(huán)境源和損傷源控制，通過盲源分離（Blind Source Separation，BSS）將二者分開后，可以根據(jù)損傷源信號是否發(fā)生突變識別損傷。刁延松和任紅［17］對結構響應數(shù)據(jù)建立自回歸模型，利用因子分析（Factor Analysis）去除溫度對自回歸系數(shù)的影響。Santos等［18-19］采用改進的高斯混合模型學習結構模態(tài)頻率在不同環(huán)境條件下的概率分布規(guī)律，利用聚類后的簇群表示不同環(huán)境因素的影響。邱雷等［20］利用GMM建立導波監(jiān)測特征參數(shù)在環(huán)境因素影響下的概率模型，采用KL距離衡量基準GMM和待測GMM的遷移趨勢，從而實現(xiàn)損傷預警。

近幾年，基于協(xié)整理論的損傷預警方法得到了快速發(fā)展。Cross等［21-22］最早對結構健康監(jiān)測中的非平穩(wěn)數(shù)據(jù)進行協(xié)整研究，發(fā)現(xiàn)模態(tài)頻率存在的協(xié)整關系不會因為環(huán)境條件的改變而發(fā)生變化，但損傷的出現(xiàn)會破壞這種均衡關系，通過協(xié)整殘差可以對損傷進行定量。刁延松等［23］提出基于自回歸模型系數(shù)和協(xié)整的結構損傷預警方法，并通過海洋平臺振動臺模型驗證該方法的有效性。Huang等［24］將結構模態(tài)頻率之間的協(xié)整系數(shù)作為卡爾曼濾波的狀態(tài)向量，通過觀察協(xié)整關系的改變實現(xiàn)在線損傷預警。李秀娟等［25］將協(xié)整理論引入到壓電阻抗法中，對溫度變化影響下的阻抗譜峰值頻率進行協(xié)整，利用協(xié)整殘差的突變識別損傷。基于協(xié)整理論的方法從非平穩(wěn)時間序列分析入手，利用協(xié)整關系消除環(huán)境因素引起的數(shù)據(jù)非平穩(wěn)性，相比于傳統(tǒng)復雜的數(shù)據(jù)分析算法具有原理簡單、易于實現(xiàn)的特點。然而，實際工程中的協(xié)整關系往往具有一定的非線性［21］，采用傳統(tǒng)的線性協(xié)整理論往往無法描述真實的協(xié)整關系，因此需要發(fā)展基于非線性協(xié)整的損傷預警方法［26-27］。

本文結合高斯混合模型的聚類算法和協(xié)整理論進行非線性環(huán)境因素影響下的損傷預警研究。該方法以結構模態(tài)頻率作為協(xié)整變量，利用高斯混合模型對不同環(huán)境條件下的頻率樣本進行概率分布擬合，然后依據(jù)高斯后驗概率將這些樣本分為不同的簇群；根據(jù)分段線性化思想，如果簇群的數(shù)量（聚類數(shù)）取值合理，那么相同簇群里的樣本近似滿足線性協(xié)整關系，采用Johansen檢驗可求出各個簇群的協(xié)整方程；對于待測樣本，同樣采用高斯后驗概率確定其簇群歸屬，之后代入相應的協(xié)整方程求出協(xié)整殘差，通過X-bar控制圖［28］實現(xiàn)結構損傷預警。

1 高斯混合聚類

高斯混合模型是高斯模型的簡單擴展［29］。假設隨機向量x∈Rl×1服從多元高斯分布，其概率密度函數(shù)可由均值向量μ∈Rl×1和協(xié)方差矩陣Σ∈Rl×l決定，具體表達式為

由此可以定義GMM的概率密度函數(shù)為

式中m代表高斯分量的個數(shù)；αi代表第i個高斯分量的混合系數(shù)，滿足代表第i個高斯分量的概率密度函數(shù)，其均值向量和協(xié)方差矩陣分別為μi和Σi。

GMM聚類算法［29］采用多個高斯分布的線性組合作為數(shù)據(jù)分布的概率密度函數(shù)，通過概率模型對應的后驗概率確定簇群的劃分。如果樣本集X=[x1x2…xs]由式（2）所表示的GMM生成，引入新的隨機變量zj∈{1，2，…，m}代表生成樣本xj的高斯分量，顯然其先驗概率p(zj=i)=αi。根據(jù)貝葉斯定理，zj的后驗概率可以表示為

采用期望最大（Expectation Maximization，EM）算法求出GMM的未知參數(shù)后，可根據(jù)式（3）計算出后驗概率最大的高斯分量作為xj的歸屬。相比于傳統(tǒng)的K均值聚類，GMM聚類產(chǎn)生的簇群形狀可以是任意的橢圓，實際應用范圍更廣。

2 協(xié)整理論

協(xié)整最早用于解決非平穩(wěn)時間序列建模時引起的虛假回歸問題。其核心思想是利用多個非平穩(wěn)變量之間存在的長期均衡（協(xié)整）關系，通過線性組合的方式將非平穩(wěn)性時間序列轉化為平穩(wěn)時間序列。

如果非平穩(wěn)時間序列yt(t=1，2，…)經(jīng)過d次差分之后剛好成為平穩(wěn)時間序列，則稱yt具有d階單整性，記為yt～I(d)。顯然，平穩(wěn)時間序列的單整性為I(0)。假定隨機向量yt={y1t y2t…ylt}T的每個分量都具有1階單整性［30］，而且存在l維向量β={β1β2…βl}T使得yt的線性組合變?yōu)槠椒€(wěn)時間序列，即βTyt～I(0)，那么yt存在協(xié)整關系，其協(xié)整方程表示為

式中β為協(xié)整向量；εt為協(xié)整殘差。

目前，協(xié)整性檢驗主要分為EG（Engle-Granger）兩步法和Johanson法。EG檢驗［31］是以殘差為基礎的檢驗方法，首先采用最小二乘法對變量進行回歸建模，再通過回歸殘差的平穩(wěn)性檢驗確定協(xié)整關系是否成立。以隨機向量yt為例，選擇y1t作為因變量建立如下回歸方程

Johanson檢驗［31］是以向量誤差修正（Vector Error Correction，VEC）模型為基礎的多變量協(xié)整性檢驗方法。VEC模型不僅包含變量間的長期協(xié)整關系，而且考慮了變量短期波動對協(xié)整方程的影響。因此，Johanson法在協(xié)整向量的估計精度上要高于采用簡單回歸分析的EG法。仍以隨機向量yt建立如下VEC模型

式中Π∈Rl×l代表長期均衡矩陣，與變量協(xié)整關系相關；Ψi∈Rl×l代表短期動態(tài)矩陣，與變量瞬時波動特征相關；Δyt∈Rl×1為yt的差分項；p為滯后階次；ut為l維白噪聲向量。矩陣Π由修正矩陣A∈Rl×r（協(xié)整殘差對Δyt的修正速度）和協(xié)整向量矩陣B∈Rl×r（每一列代表一個協(xié)整向量）組成，其中r代表矩陣Π的秩（r=0說明yt不存在協(xié)整關系）。因為yt～I(1)，則Δyt，Δyt-i～I(0)，式（6）中的所有變量都具有平穩(wěn)性，采用極大似然法估計協(xié)整向量矩陣B，然后選擇第一列協(xié)整向量（對應特征值最大、平穩(wěn)性最強［31］）建立協(xié)整方程。

3 結構損傷預警方法

結構各階模態(tài)頻率受環(huán)境因素影響存在協(xié)整關系［30］，損傷的出現(xiàn)將破壞這種關系而導致協(xié)整殘差出現(xiàn)突變，通過比較損傷前后的協(xié)整殘差可以實現(xiàn)損傷預警。但是，環(huán)境因素影響往往存在一定的非線性，協(xié)整變量間難以保證較好的線性協(xié)整關系，導致基于線性協(xié)整理論的損傷預警方法精度不高。為解決上述問題，本文根據(jù)分段線性化思想，利用GMM聚類將非線性協(xié)整關系轉化為多個線性協(xié)整關系，彌補協(xié)整理論無法處理非線性環(huán)境因素影響的問題。

3.1 幾何解釋

為說明GMM聚類結合線性協(xié)整理論具有處理非線性環(huán)境因素影響的能力，以模態(tài)頻率f1和f2組成的二維樣本為例。如圖1所示，由于環(huán)境因素的影響，f1和f2存在明顯的雙線性關系［32］，線性協(xié)整理論無法準確描述二者的非線性協(xié)整關系。本文采用GMM聚類將這些頻率樣本分割成簇群1和簇群2。根據(jù)分段線性化思想，假定聚類后每個簇群都不存在非線性特征，采用Johanson法建立兩個線性協(xié)整方程近似原來的非線性協(xié)整關系。

圖1 非線性環(huán)境因素影響下模態(tài)頻率f1和f2的散點分布圖Fig.1 Scatter plot of modal frequencies f1 and f2 under nonlinear environmental effects

對于損傷樣本點A0，計算其后驗概率可確定其簇群歸屬：如果A0屬于簇群1的概率較大，那么利用協(xié)整方程1求出協(xié)整殘差（線段A0A1的長度）作為損傷指標；反之，如果A0屬于簇群2的概率較大，那么根據(jù)協(xié)整方程2求得協(xié)整殘差（線段A0A2的長度）作為損傷指標。

3.2 損傷預警流程

土木工程結構可能產(chǎn)生的損傷千差萬別，而且在服役期間的大部分時段都處于正常狀態(tài)，缺少或很難建立起完整的損傷樣本庫，無法通過樣本匹配的方式確定未知樣本是否屬于某類損傷樣本。最常用的方法是將損傷樣本看作離群點，采用X-bar控制圖［28］實現(xiàn)損傷預警。針對正態(tài)分布的樣本總體，X-bar控制圖根據(jù)小概率事件原理確定警戒線，對異常事件進行實時監(jiān)測和報警。

假定結構無損傷（參考）樣本的協(xié)整殘差滿足正態(tài)分布，根據(jù)99%置信度確定控制上限（Upper Control Limit，UCL）和 控 制 下 限（Lower Control Limit，LCL）分別為

式中CL和σ分別代表協(xié)整殘差的均值和標準差。當待測樣本的協(xié)整殘差落在上下控制限范圍之外，則說明結構出現(xiàn)損傷。

因為GMM聚類得到的不同簇群可能存在不同的離散程度（即不同簇群的樣本標準差σ存在差異），如圖1中簇群2的樣本離散性大于簇群1，導致根據(jù)式（7）確定的控制上下限也會出現(xiàn)不同：采用簇群2對應的X-bar控制圖進行損傷預警，由于其控制上下限之間的區(qū)域較寬，某些歸屬于簇群1的損傷樣本可能會被誤認為無損樣本；反之，采用簇群1對應的X-bar控制圖會導致某些歸屬于簇群2的健康樣本被誤判為損傷樣本。因此，需要消除不同簇群的樣本離散性對X-bar控制圖的影響，最簡單的方法是對歸屬于簇群i的協(xié)整殘差樣本進行標準化處理，處理后的協(xié)整殘差可表示為

式中σ()i為簇群i所有參考樣本的協(xié)整殘差標準差。理論上，標準化后所有參考樣本的協(xié)整殘差都滿足標準正態(tài)分布N（0，1），可采用統(tǒng)一的X-bar控制圖對不同簇群的待測樣本進行監(jiān)測。

綜上，基于GMM聚類和協(xié)整理論（簡寫為GMM-CI）的損傷預警方法分為以下4步。

Step 1：模態(tài)參數(shù)識別。采集無損傷結構在不同環(huán)境條件下的振動響應數(shù)據(jù)，并采用運行模態(tài)分析技術［33］（如隨機子空間法、頻域分解法等）識別結構模態(tài)頻率作為協(xié)整變量。

Step 2：GMM聚類。將Step 1得到的頻率樣本作為參考樣本，利用GMM聚類對其進行聚類分析。其中，聚類數(shù)m，可根據(jù)貝葉斯信息準則［34］（Bayesian Information Criterion，BIC）確定，其表達式為

式中L（m）為似然函數(shù)；n(m)為待估計參數(shù)個數(shù)；s為參考樣本數(shù)量。隨著聚類數(shù)k的增加，BIC值會先減少后增大，選擇BIC最小值對應的m值作為最佳聚類數(shù)。

Step 3：估計協(xié)整方程。對聚類后相同簇群的參考樣本進行Johansen檢驗，選擇最大特征值對應的協(xié)整向量建立如式（4）所示的協(xié)整方程，同時獲得該簇群的協(xié)整殘差樣本。

Step 4：建立X-bar控制圖。計算每個簇群協(xié)整殘差樣本的標準差σ()i，并按照式（8）對樣本進行標準化處理，之后根據(jù)式（7）確定X-bar控制圖的控制上下限。

Step 5：損傷預警。對于結構未知狀態(tài)的振動響應，首先按照Step 1得到其頻率樣本，然后通過式（3）計算GMM后驗概率確定簇群歸屬，之后利用對應的協(xié)整方程求出協(xié)整殘差，最后根據(jù)式（8）對協(xié)整殘差進行標準化處理，如果標準化后的協(xié)整殘差沒有落在X-bar控制圖上下限的范圍內(nèi)，則說明結構出現(xiàn)損傷。

4 試驗數(shù)據(jù)驗證

為驗證上述方法的有效性，將其應用于瑞士Z24橋的現(xiàn)場試驗數(shù)據(jù)。該橋為預應力混凝土箱梁橋，主跨30 m，兩邊跨均為14 m，如圖2（a）所示。為研究環(huán)境因素以及損傷對結構動力參數(shù)的影響，從1997年11月11日 到1998年9月11日，Z24橋 安 裝的健康監(jiān)測系統(tǒng)采集了結構加速度響應數(shù)據(jù)和包括環(huán)境溫濕度、風速、風向等在內(nèi)的環(huán)境數(shù)據(jù)。布置的16個加速度傳感器在監(jiān)測過程中只有8個（編號3，5，6，7，10，12，14，16）保持正常工作，如圖2（b）和（c）所示。結構漸進破壞試驗持續(xù)了一個月（1998年8月9日到1998年9月9日），設置的損傷［32］依次是橋墩沉降、混凝土脫落、墩臺滑坡、混凝土鉸接失效、錨頭失效和鋼筋破裂。

圖2 瑞士Z24橋及其加速度傳感器分布［32］Fig.2 Switzerland Z24 bridge and the layout of acceleration sensors［32］

Peeters和De Roeck［6］采 用 隨 機 子 空 間 法 識 別Z24橋的前4階模態(tài)頻率，其變化如圖3所示。因為環(huán)境因素的影響，即使結構沒有出現(xiàn)損傷（樣本點1到3470），其模態(tài)頻率仍出現(xiàn)不同程度的波動，而且高階模態(tài)頻率的波動性明顯大于低階模態(tài)。漸進破壞試驗（樣本點3471到3932）導致結構各階模態(tài)頻率均呈現(xiàn)下降的趨勢。

圖3 瑞士Z24橋前4階模態(tài)頻率的變化趨勢Fig.3 Variations of the first four modal frequencies of the Z24 Bridge

Z24橋的模態(tài)頻率隨著橋梁路面瀝青層溫度的升高而降低，同時存在明顯的雙線性關系（如圖4所示），其原因［32］可解釋為：當溫度大于0℃時，路面瀝青的彈性模量改變對結構剛度的影響起主導作用，溫度越高彈性模量越低，導致結構模態(tài)頻率隨著溫度升高而降低；當溫度小于0℃時，結冰引起的結構邊界條件改變起主導作用，低溫增強了邊界約束，同時也增大結構剛度，出現(xiàn)結構模態(tài)頻率隨著溫度降低而增大的現(xiàn)象。從協(xié)整理論的角度看，以模態(tài)頻率作為協(xié)整變量至少存在兩種均衡關系，是一種非線性協(xié)整關系。

圖4 結構模態(tài)頻率和溫度的雙線性關系Fig.4 Bi-linear relationship between structural modal frequency and temperature

4.1 損傷預警結果

以f1和f2作為協(xié)整變量，將圖3中編號1到3000的樣本作為參考（訓練）樣本，編號3001到3932的樣本作為待測樣本（其中編號3001到3470的樣本是無損傷樣本，編號3471到3932的樣本是損傷樣本），采用本文所提基于高斯混合聚類和協(xié)整理論的損傷預警方法預報結構損傷。

首先采用BIC準則確定高斯混合聚類的最佳聚類數(shù)。如圖5所示，聚類數(shù)為3時的BIC值達到最小。因此，采用包含3個高斯分量的GMM對參考樣本進行聚類，其結果如圖6所示。根據(jù)各簇群樣本對應的環(huán)境溫度數(shù)據(jù)（見圖7）可知，簇群1（對應的環(huán)境溫度大于0℃）和簇群3（對應的環(huán)境溫度小于0℃）代表上述兩種由溫度主導的均衡機制；簇群2（既有大于0℃的樣本，也有小于0℃的樣本）可能是其他環(huán)境因素（風速或車載等）主導的均衡機制。

圖5 BIC準則確定GMM聚類數(shù)Fig.5 BIC criterion to determine the cluster number

圖6 參考樣本GMM聚類結果Fig.6 GMM clustering results of reference samples

圖7 各簇群模態(tài)頻率與溫度對應關系Fig.7 Relationship between modal frequencies and temperature in each cluster

對各簇群的參考樣本進行ADF檢驗，設定置信水平為5%置信水平，則p＜5%（或者t-統(tǒng)計量小于-1.941）代表變量是平穩(wěn)時間序列。從表1可出，模態(tài)頻率f1和f2是非平穩(wěn)時間序列，但其一階差分是平穩(wěn)時間序列，即模態(tài)頻率f1和f2都具有一階單整性。

表1 模態(tài)頻率f1和f2及其一階差分的ADF檢驗結果（5%置信水平）Tab.1 ADF test for modal frequencies f1，f2 and their first differences（5% confidence level）

確定模態(tài)頻率f1和f2都滿足一階單整條件后，采用Johanson檢驗計算各簇群對應的協(xié)整方程，其表達式為

利用上式可得不同簇群參考樣本的協(xié)整殘差，之后根據(jù)式（8）對協(xié)整殘差進行標準化處理，最后確定X-bar控制圖的控制上下限。

對于待測樣本，根據(jù)式（3）計算其屬于不同簇群的概率，選擇最大概率值對應的簇群作為最終歸屬，之后將樣本代入式（10）中對應簇群的協(xié)整方程，求得協(xié)整殘差并按照式（8）對其進行標準化處理，最后通過X-bar控制圖進行損傷預警，結果如圖8所示。從圖中可以看出：參考樣本的協(xié)整殘差不存在明顯的變化趨勢（樣本點1500附近的異常突變與圖6中畫圈部分樣本的模態(tài)識別偏差有關），說明采用本文所提方法能夠有效消除非線性環(huán)境因素的影響；待測損傷樣本的協(xié)整殘差隨損傷的累積有變大的趨勢，說明協(xié)整殘差能體現(xiàn)損傷的相對嚴重程度。需要指出的是，待測樣本中有4.89%的無損傷樣本被誤判為損傷樣本，2.38%的損傷樣本被誤判為無損傷樣本。為進一步分析造成上述誤差的原因，圖9給出了歸屬于簇群1的參考樣本以及所有待測樣本的分布情況（待測樣本的聚類結果均為簇群1）。從圖中可以看出，出現(xiàn)誤判的待測樣本主要位于警戒線2附近：對于待測無損傷樣本，理論上99%置信度確定的X-bar控制圖會有1%的樣本均勻分布在兩條警戒線之外，而實際位于警戒線2以下的樣本有19個，位于警戒線1以上的樣本只有4個，其原因最有可能是局部樣本的頻率識別誤差過大，導致位于警戒線2以下的誤判樣本遠超過位于警戒線1以上的樣本；至于損傷樣本，由于漸進破壞試驗初期損傷引起的協(xié)整殘差變化不明顯，再加上測量噪聲引起的樣本不確定性，導致個別損傷樣本被誤認為無損傷樣本。

圖8 基于GMM-CI方法的損傷預警結果（以f1和f2作為協(xié)整變量）Fig.8 Damage alert results based on GMM-CI method（using f1 and f2 as cointegration variables）

圖9 參考樣本與待測樣本的散點分布Fig.9 Scatter distribution of reference sample and samples to be tested

4.2 不同協(xié)整方法比較結果

為比較不同協(xié)整方法的結果，圖10給出了傳統(tǒng)基于線性協(xié)整理論的損傷預警結果。由于線性協(xié)整模型無法完全消除非線性環(huán)境因素的影響，損傷樣本的誤判率高達41.77%，在識別精度上遠低于本文所提的基于GMM聚類相結合的協(xié)整理論的損傷預警方法。

圖10 基于線性協(xié)整理論的損傷預警結果Fig.10 Damage alert using linear co-integration theory

圖11為Shi等［31］提出的基于區(qū)制轉移的非線性協(xié)整方法的損傷預警結果。該方法利用環(huán)境溫度信息識別出模態(tài)頻率之間存在兩種協(xié)整關系（以0.98℃為界），通過待測樣本對應的溫度信息劃分其協(xié)整關系歸屬，之后根據(jù)對應的協(xié)整方程求得協(xié)整殘差作為損傷指標。由于考慮了溫度信息，基于區(qū)制轉移的非線性協(xié)整方法的誤判率（無損傷樣本為2.77%，損傷樣本為1.51%）均低于本文所提方法。但是，考慮環(huán)境因素的監(jiān)測難度以及成本控制（比如邊界條件、車載等），本文所提方法無需借助環(huán)境信息，更易于實際應用。

圖11 基于區(qū)制轉移的非線性協(xié)整方法得到的協(xié)整殘差（利用f1-f4和溫度數(shù)據(jù)）Fig.11 Cointegration residuals based on regime-switching cointegration approach（using f1-f4 and temperature data）

4.3 聚類個數(shù)對損傷預警結果的影響

BIC準則作為理想的模型選擇方法，在實際應用中易受樣本分布和測量噪聲的影響，很有可能產(chǎn)生錯誤的結果，有必要研究不同聚類數(shù)對損傷預警精度的影響，并提出合理的聚類數(shù)建議。如圖12所示，當聚類數(shù)取2時，待測損傷樣本的誤判率較大，出現(xiàn)欠擬合現(xiàn)象；當聚類數(shù)取3，4和5時，損傷預警的誤判率均處于5%附近。顯然，欠擬合會大大降低本文所提方法的損傷預警能力，必須設法避免聚類數(shù)取值過小導致?lián)p傷預警失效。

圖12 不同聚類數(shù)對損傷預警精度的影響Fig.12 Influences of different cluster numbers on damage alert accuracy

考慮GMM聚類的目的是將非線性協(xié)整關系轉化為多個線性協(xié)整關系的組合，因此聚類數(shù)的選擇理論上應保證聚類后各個簇群都不存在非線性特征。圖13（a）中，簇群2由于聚類數(shù)不足仍具有非線性特征，采用線性協(xié)整理論進行建模必然會有誤差。根據(jù)分段線性化思想，簇群數(shù)越多（分段數(shù)越多）越能近似非線性特征，比如圖13（b）簇群1，4和5相比于圖13（a）簇群1有更多細節(jié)上的線性特征，因此聚類數(shù)取4或5仍能較為準確地預警損傷。但是，聚類數(shù)取值太大不僅會造成過擬合現(xiàn)象，而且嚴重影響GMM聚類的計算效率和精度。所以本文建議在BIC準則的基礎上，借助交叉驗證［35］方法產(chǎn)生多個備選聚類數(shù)，選擇較大者作為最終聚類數(shù)。

圖13 不同聚類數(shù)對GMM聚類結果的影響Fig.13 GMM clustering results with different cluster numbers

4.4 多協(xié)整變量的損傷預警

高維數(shù)據(jù)一般存在多種協(xié)整關系，不同的協(xié)整關系對損傷預警精度的影響很大。實際應用中可以采 用 主成分分 析［36］（Principle Component Analysis，PCA）將高維數(shù)據(jù)進行降維處理，通常以第一主成分和第二主成分作為新的協(xié)整變量，利用二維數(shù)據(jù)協(xié)整關系的唯一性特點避開高維數(shù)據(jù)協(xié)整關系的選擇難題。除此之外，降維不僅能夠簡化GMM聚類的復雜度，而且二維空間便于人工判斷BIC準則確定的聚類數(shù)以及GMM聚類結果是否合理。

以 模 態(tài) 頻 率f1，f2和f3為 例，同 樣 選 擇 圖3前3000個數(shù)據(jù)作為參考樣本，采用PCA對數(shù)據(jù)進行降維處理，將得到的第一主成分P1和第二主成分P2作為新的協(xié)整變量，然后采用GMM聚類對降維后的參考樣本進行聚類，其結果如圖14所示。相比于圖6，其簇群分布顯得更加緊湊，但依然能夠區(qū)分出3種協(xié)整關系。對于待測樣本，同樣需要經(jīng)過PCA降維，再確定其簇群歸屬。

圖14 PCA降維后參考樣本的GMM聚類結果Fig.14 GMM clustering results of reference samples after PCA

表2為各聚類參考樣本的ADF檢驗結果。從表中可以看出，主成分P1和P2是非平穩(wěn)時間序列，但其一階差分是平穩(wěn)時間序列，說明PCA分析得到的P1和P2都具有一階單整性。

表2 主成分P1和P2及其一階差分的ADF檢驗結果（5%置信水平）Tab.2 ADF test for principle components P1，P2 and their first differences（5% confidence level）

確定主成分P1和P2都滿足一階單整條件后，采用Johanson檢驗計算出各簇群對應的協(xié)整方程，其表達式為

之后的損傷預警參考4.1節(jié)。

圖15和16分別給出了未采用PCA和采用PCA后本文所提方法的損傷預警結果。前者以最小特征值對應的協(xié)整向量建立各簇群的協(xié)整方程，無損傷樣本和損傷樣本的誤判率分別為3.40%和2.60%；后者利用PCA降維技術避開協(xié)整向量的選擇問題，無損傷樣本和損傷樣本的誤判率分別為2.98%和3.03%。二者誤判率相差不大，說明采用PCA對多維協(xié)整變量進行降維處理對損傷預警精度的影響不大。

圖15 基于GMM-CI方法的損傷預警結果（最小特征值對應的協(xié)整向量）Fig.15 Damage alert results based on GMM-CI method（using the cointegration vector corresponding to the smallest eigenvalue）

圖16 基于PCA-GMM-CI方法的損傷預警結果Fig.16 Damage alert results based on PCA-GMM-CI method

最后，圖17給出了基于線性協(xié)整方法的損傷預警結果。由于未采用GMM聚類對非線性協(xié)整關系進行線性化處理，其損傷樣本的誤判率仍有42.64%。需要說明的是，無論是基于GMM-CI的方法，還是基于線性協(xié)整的方法，都存在協(xié)整向量的選擇問題。如果選擇最大特征值對應的協(xié)整方程進行損傷預警，兩種方法的損傷樣本誤判率均高達90%以上。其原因可解釋為：雖然最大特征值對應的協(xié)整方程平穩(wěn)性最強，但也說明這種協(xié)整狀態(tài)很難被損傷打破（或者說對損傷不敏感），比如基于線性協(xié)整的方法中對應于最大特征值的協(xié)整方程為

圖17 基于線性協(xié)整方法的損傷預警結果（最小特征值對應的協(xié)整向量）Fig.17 Damage alert results based on linear cointegration method（using the cointegration vector corresponding to the smallest eigenvalue）

根據(jù)圖3可知，f2對損傷最為敏感，但式（12）中其權重只有0.73，導致?lián)p傷引起的協(xié)整殘差εt變化過小而被噪聲覆蓋。

5 結論

為解決傳統(tǒng)線性協(xié)整理論無法消除非線性環(huán)境因素影響的問題，本文將高斯混合聚類引入到基于協(xié)整的損傷預警方法。該方法以結構模態(tài)頻率的協(xié)整關系為基礎，采用高斯混合聚類對存在非線性協(xié)整關系的協(xié)整變量進行分段線性化處理，聚類后相同簇群的樣本便可采用線性協(xié)整理論進行建模。對于待測樣本，先通過計算GMM后驗概率確定其簇群歸屬，再根據(jù)對應的協(xié)整方程求得協(xié)整殘差，最后利用X-bar控制圖監(jiān)測協(xié)整殘差是否出現(xiàn)異常進行損失預警。

通過Z24橋的長期監(jiān)測數(shù)據(jù)和破壞試驗對上述方法進行驗證，得到以下結論：（1）采用高斯混合聚類將協(xié)整變量間的非線性關系轉化為多個線性關系的集合，從而將基于協(xié)整理論的損傷預警方法擴展到非線性協(xié)整理論；（2）利用主成分分析對協(xié)整變量進行降維處理可以在保證損傷預警精度的前提下，避開高維數(shù)據(jù)協(xié)整關系的選擇問題，同時簡化高斯混合模型的復雜度；（3）本文所提方法不需要測量環(huán)境信息，能夠有效消除實際工程中非線性環(huán)境因素對損傷預警的影響，降低損傷預警的誤判率。

GMM聚類采用EM算法估計模型參數(shù)，聚類結果易受EM算法的初值影響，導致?lián)p傷預警結果產(chǎn)生顯著變化，如何提高GMM聚類的精度和穩(wěn)定性是未來重點研究方向。

致謝感謝比利時魯汶大學的Guido De Roeck教授提供Z24橋的模態(tài)頻率和環(huán)境溫度數(shù)據(jù)。