基于MSVAR 模型的非線性損傷識別方法

鄭 泓，段忠東，歐進萍

(哈爾濱工業大學(深圳)土木與環境工程學院，深圳 518055)

土木工程結構(如大跨橋梁、超高層建筑等)在服役過程中難免會出現損傷，嚴重時甚至會導致結構失效而引發災難性事故，因此有必要對結構的損傷狀態進行監測[1-4]。傳統的無損檢測技術需要預知損傷位置，對檢測設備難以到達的位置存在盲點，而且過于依賴巡檢人員的既有經驗和主觀判斷[3]。為克服傳統方法的不足，基于振動的損傷識別方法開始受到學者們的關注[4-5]。這類方法通常需要借助線性模型對結構振動響應進行建模，提取與損傷相關的模型參數(如模態頻率等)作為損傷敏感特征，通過比較損傷前后這些特征的變化實現損傷識別。

實際工程中出現的損傷往往具有一定的非線性特征[6-7]，比如裂縫的開合效應導致結構響應出現非平穩特征[8]；結構在地震作用下的塑性耗能使得恢復力曲線形成近似封閉的滯回環[9]。采用線性模型不僅無法獲取損傷的真實信息，而且極易受到其他非線性因素的干擾，導致損傷預警的精度和可信度不高。以模態頻率為例，環境溫度變化對模態頻率的影響往往大于損傷的影響[10-12]。因此，有必要引入非線性分析的方法，通過準確獲取損傷的非線性信息提高損傷識別方法的工程適用性。

非線性損傷識別主要分為非線性指標方法和非線性系統識別方法兩大類[6]。其中，基于非線性指標的方法認為損傷會導致結構響應出現非線性特征，借助能夠表現數據非線性變化的指標實現損傷預警。該方法具有建模方便，易于實現在線監測的優點，已成為非線性損傷識別的研究熱點。Robertson 和Farrar 等[13]認為非線性損傷的出現會導致結構響應產生突變，利用Holder 指數跟蹤結構響應信號的不連續性實現損傷監測。Chen和Yu[14]結合自回歸滑動平均模型(ARMA)和廣義自回歸條件異方差模型(GARCH)對結構加速度響應進行建模，定義模型殘差的標準差作為損傷敏感特征，通過試驗表明采用非平穩時間序列模型能夠增強損傷識別的準確性。Yu 和Zhu[15]對結構振動響應建立自回歸模型(AR)，利用模型殘差的高階矩識別非線性損傷。馬家欣和許飛云等[16]將帶有外部輸入的非線性自回歸模型(GNARX)應用于鋼板的損傷識別，發現其識別效果優于線性時間序列模型。Li 和Sun 等[17]通過計算結構自由響應信號的互相干矩陣識別結構非線性損傷，克服了傳統基于相干函數的非線性損傷識別方法需要測量激勵信號的缺點。Chen 和Yu 等[18]為解決ARMA 模型無法處理非線性損傷識別的問題，構造向量空間余弦相似度作為損傷敏感特征，然后利用K 均值聚類和貝葉斯判別實現非線性損傷的快速預警。

在非線性系統識別領域中，損傷模型用于描述不同類型損傷對結構物理特性(主要是剛度和阻尼)的影響，比如土木工程結構中較為常見的損傷裂縫用雙線性剛度模型表示。上述基于非線性指標的方法把關注點放在“是否出現非線性損傷”，忽略了損傷非線性特征的差異性，無法識別損傷類型。基于非線性系統識別的方法采用不同損傷模型對結構響應進行建模，以最小建模誤差作為確定損傷類型的標準。該方法在建模上需要消耗大量時間[19]，不適合結構在線監測，而且識別損傷類型之前仍然需要判斷結構是否出現非線性損傷。因此，本文結合非線性指標方法和非線性系統識別方法的優點，提出基于馬爾科夫區制轉移向量自回歸模型(MSVAR)的非線性損傷識別方法，并通過數值算例和模型試驗驗證本文所提方法能夠有效識別裂縫損傷。

1 MSVAR 模型理論

MSVAR 模型屬于非線性時間序列模型，最早用于挖掘經濟行為的內在機制，比如國民生產總值受政策影響而產生的增長期和衰減期[20]。

對于滿足零均值但非平穩的l維時間序列Yt∈Rl×1(t=1,2,···,T)，MSVAR 模型描述這樣一個數據過程[21]：Yt受m個機制影響，相同機制下的數據用同一個p階向量自回歸模型(VAR)描述；數據在不同機制之間的轉換采用馬爾科夫鏈表示。上述MSVAR 模型簡寫成MS(m)VAR(p)，用公式表示為：

MSVAR 模型采用隱狀態平滑概率表征數據過程處于某一種隱狀態的可能性大小，在實際應用中可通過隱狀態平滑概率的改變判斷研究對象在何時出現狀態改變。在計算平滑概率前，需要采用極大似然法估計MSVAR 模型參數[21]。假設式(1)中Yt取自隱狀態i的條件概率滿足多維正態分布，即：

式中：符號 (÷) 代表兩個向量對應元素相除；ξt|T的第i個元素代表第i個隱狀態在t時刻出現的概率。

2 非線性損傷預警原理

非線性損傷指標法假設損傷會導致結構振動響應出現非線性特征，即結構從線性系統轉變成非線性系統的預示著損傷的出現。考慮到環境激勵下線性系統的振動響應可用VAR 模型表示[22-24]，所以采用MSVAR 模型對響應數據進行建模會出現兩種情況：第一種是MSVAR 模型某一隱狀態的平滑概率始終為100%，表明數據過程只受一種機制作用，MSVAR 模型退化為線性的VAR 模型，結構不存在損傷；第二種是MSVAR 模型有多個隱狀態的平滑概率大于零，說明采用非線性的MSVAR 模型才能準確描述結構振動特征，即結構出現損傷。

采用四階龍格庫塔法計算振子在白噪聲激勵下的響應，以100 Hz 采樣頻率采集其位移數據進行分析。如圖1 所示，該單自由度振子由小幅簡諧振動過度到幅值較大的簡諧振動。其中，過度階段(2 s～3 s)響應幅值逐漸增大的原因是 |x|＜1導致系統阻尼為負。

圖1 線性系統到非線性系統的位移響應時程Fig. 1 Time history of displacement response from linear system to nonlinear system

對生成的400 個位移響應數據建立MS(3)VAR(2)模型，其隱狀態平滑概率時程如圖2 所示。前2 s 隱狀態2 和隱狀態3 的平滑概率均為0，MSVAR 模型退化為隱狀態1 對應的VAR 模型，說明結構響應可以采用線性模型進行建模；2 s后隱狀態1 的平滑概率為0，隱狀態2 和隱狀態3的平滑概率都存在大于零的時刻點，說明采用簡單的線性VAR 模型已經不能準確描述結構真實振動狀態，結構由線性系統轉為非線性系統。

圖2 線性系統到非線性系統的平滑概率變化Fig. 2 Changes of smooth probabilities from linear system to nonlinear system

為實現損傷預警，本文引入信息熵作為非線性損傷的預警指標，其表達式[27]為：

3 裂縫類型損傷的識別方法

判斷結構出現損傷后需要對損傷類型進行識別，并進一步確定損傷位置和程度。因為不同類型的損傷對應不同的數學模型[6]，本節根據MSVAR 模型與雙線性剛度系統的等價關系，闡述識別裂縫類型損傷的思想。

3.1 雙線性剛度系統的MSVAR 模型

雙線性剛度系統最極端的例子是碰撞振子，該系統對應于一個球對著一面剛性墻壁的碰撞。在土木工程領域中，與碰撞振子類似的是帶有裂縫損傷的結構[28]，以單自由度系統為例，其運動方程為：

3.2 損傷類型識別

裂縫的張開和閉合導致結構在兩種不同的振動狀態來回切換，對應的MSVAR 模型也會出現兩個隱狀態交替激活的現象，即平滑概率時程在0 和100%兩個值之間切換。

以3 自由度矩形截面懸臂梁為例展示裂縫損傷對隱狀態平滑概率影響的特殊性，其振動微分方程為：

根據截面慣性矩定義可知，式(15)表明裂縫張開時裂縫位置懸臂梁的有效橫截面積約為裂縫閉合時的85%。

同樣，采用四階龍格庫塔法計算上述懸臂梁在白噪聲激勵下的振動響應，以50 Hz 采樣頻率收集200 個位移數據建立MS(3)VAR(2)模型，各隱狀態的平滑概率時程如圖3 所示。從圖3 可以看出，MSVAR 模型只有隱狀態1 和隱狀態2 被激活，與裂縫張開和閉合兩種結構振動狀態相符；隱狀態平滑概率大部分在0 和100%兩個值之間變化，滿足裂縫損傷特有的變化模式。

圖3 懸臂梁平滑概率時程Fig. 3 Time histories of smooth probability for cantilever beam

為計算方便，這里選擇平滑概率由0 突變到100%作為裂縫的損傷模式，通過統計其出現次數N0→1判斷損傷類型是否屬于裂縫。N0→1值越大表明出現裂縫損傷的可能性也越大，比如范德波爾系統(圖2)和3 自由度懸臂梁(圖3)的N0→1值分別為1 和10。

3.3 損傷定位

圖4 懸臂梁損傷定位結果Fig. 4 Damage location result for cantilever beam

3.4 損傷定量

裂縫損傷的嚴重程度一般用裂縫寬度表示。由式(15)可知，裂縫寬度可用裂縫開合瞬間(設為ts)損傷單元處的層間位移 |x2(ts)-x3(ts)|表示。根據3.2 節結論，這里的裂縫開合瞬間可代替為裂縫損傷模式出現的時間點tsi(i=1,2,···,N0→1)。因此，假設裂縫位于自由度i和j之間，其寬度Wc可表示為：

圖5 給出3.2 節3 自由度系統的損傷定量結果。從圖5 可以看出，由于存在樣本不確定性，大部分裂縫寬度的估計值與真實值1 mm 偏差較大；按照式(17)對樣本取平均可消除樣本偏差，得到較為準確的裂縫寬度估值。

圖5 懸臂梁損傷定量結果Fig. 5 Quantitative results of crack damage for cantilever beam

3.5 損傷識別流程

根據上述分析，基于MSVAR 模型的非線性損傷識別方法可分為兩部分：第一部分是對結構響應建立MSVAR 模型，通過隱狀態平滑概率構造信息熵判斷結構是否出現損傷；第二部分是根據MSVAR 模型能夠準確描述結構出現裂縫后的振動狀態，進一步對損傷類型、損傷位置和損傷程度進行識別。損傷識別流程如圖6 所示，具體可分為以下5 個步驟：

圖6 非線性損傷識別流程Fig. 6 Nonlinear damage identification flowchart

步驟1：建立MSVAR 模型。收集結構振動數據并采用擬牛頓算法[29]估計MSVAR 模型參數。建模時需要根據均方根誤差準則(RMSE)選擇合適的VAR 模型階數[23]，然后采用貝葉斯信息準則(BIC)確定隱狀態個數[21]。

步驟2：非線性損傷預警。計算式(8)隱狀態平滑概率時程，之后根據式(11)求出信息熵Ent作為損傷預警指標。如果Ent大于零，則認為結構存在損傷。

步驟3：損傷類型識別。統計隱狀態平滑概率從0 突變到100%的次數N0→1判斷出現的損傷是否屬于裂縫。N0→1值越大說明裂縫損傷的可能性越大。

步驟4：裂縫類型損傷定位。根據式(16)計算損傷定位向量dlv，以dlv最大值元素對應的自由度作為裂縫位置。

步驟5：估計裂縫寬度。根據步驟3 得到的N0→1個平滑概率突變時刻點以及步驟4 的定位結果，按照式(17)估計裂縫寬度。

4 試驗驗證

以美國洛斯阿拉莫斯實驗室三層框架結構的試驗數據[30]驗證本文所提方法的有效性。如圖7所示，該結構由鋁板( 3 0.5 cm×30.5 cm×2.5 cm)和鋁柱( 1 7.7 cm×2.5 cm×0.6 cm)通過螺栓連接而成。在結構第三層固定有懸臂柱和阻尼器，通過調節二者之間的縫隙寬度模擬裂縫損傷大小。受單向軌道約束，結構只能沿激勵方向振動，白噪聲激勵的帶寬為20 Hz～150 Hz。加速度傳感器安裝在每層鋁板中間，采樣頻率為320 Hz。

圖7 實驗室三層框架結構[30]Fig. 7 Laboratory three-layer frame structure

試驗總共設有17 個工況，如表1 所示。每個工況收集10 個加速度數據樣本，每個樣本均有320 個數據。為模擬環境噪聲影響，采集的加速度數據混有高斯白噪聲，數據信噪比為20 dB。

表1 三層框架結構的損傷工況[30]Table 1 Damage cases of three-layer frame structure

4.1 非線性損傷預警

對采集的加速度數據建立MSVAR 模型需要確定VAR 模型階數和隱狀態個數。首先采用RMSE準則選擇合適的VAR 模型階數。圖8(a)給出4 個典型工況的RMSE 值隨VAR 模型階數變化的曲線，根據曲線拐點可判斷VAR 模型階數應取為2。確定VAR 模型階數之后，需要根據BIC 準則選擇隱狀態個數。如圖8(b)所示，工況14 的BIC值在隱狀態個數為3 時達到最小，而其他工況的最佳隱狀態個數為2。為避免欠擬合現象，這里選擇MS(3)VAR(2)模型對結構加速度響應數據進行擬合。

圖8 MSVAR 模型階數確定結果Fig. 8 MSVAR model order determination results

圖9 對比了基于MS(3)VAR(2)模型和基于MS(2)VAR(2)模型的損傷預警結果。由于數據噪聲和建模誤差等因素的影響，無損傷樣本的信息熵一般不為零，需要根據參考工況(工況1)的信息熵設置損傷閾值，信息熵小于損傷閾值說明結構沒有損傷。從圖9 可以看出，兩種MSVAR 模型均能準確預警損傷工況11～工況14 和工況17；但損傷工況10、工況15 和工況16 都被誤認為健康工況，其原因是過大的縫隙寬度使得結構振動時阻尼器與懸臂柱幾乎不發生碰撞，導致計算的信息熵很難與無損傷樣本區分開。

隱狀態個數過少的缺點主要體現在工況14 的信息熵整體偏小(圖9(b)藍色橢圓位置)，甚至有個別樣本的信息熵因小于損傷閾值而被誤認為無損傷樣本。考慮到工況14 屬于小損傷工況(縫隙寬度只有0.05 mm)，隱狀態個數取值過小可能會降低本文所提方法對損傷的靈敏度。

圖9 三層框架結構損傷預警結果Fig. 9 Damage warning results for frame structure

4.2 裂縫類型損傷識別

確定結構存在損傷后，需要根據MSVAR 模型平滑概率是否具有從0 突變到100%的變化模式進一步判斷出現的損傷是否屬于裂縫類型損傷。

以工況11～工況14 為例，這些損傷工況對應的MSVAR 模型隱狀態平滑概率時程曲線如圖10所示。從圖10 可以看出，4 個損傷工況的MSVAR模型都有2 個隱狀態被激活，而且隱狀態之間的變化模式與裂縫損傷相同，存在平滑概率從0 突變到100%的特殊損傷模式。出現這種現象的原因是懸臂柱與阻尼器的碰撞使得結構第三層剛度出現雙線性特征，繼而導致結構在兩種振動狀態來回切換。

圖10 不同損傷工況平滑概率曲線的對比結果Fig. 10 Comparisons of smooth probability curves under different damage cases

考慮建模誤差和噪聲干擾，裂縫的損傷模式放寬為隱狀態平滑概率由小于5%突變到大于95%，其發生次數(仍記為N0→1)的統計結果如圖11所示。因為工況11～工況14 的縫隙寬度由寬變窄，而減少縫隙寬度意味著結構振動時阻尼器與懸臂柱的碰撞機會增多(從圖10 中也能看出隱狀態1 和隱狀態2 之間的切換頻次與縫隙寬度成反比)，所以工況14 的N0→1值最大，最有可能是裂縫類型損傷。需要特別指出的是，工況14 雖然屬于小損傷工況，但其信息熵(圖9(a))和N0→1值(圖11)都遠大于其他損傷工況，說明本文所提方法對早期裂縫較為敏感。

圖11 裂縫損傷的隱狀態轉移模式出現次數統計Fig. 11 Statistical results of occurrence number for crack damage hidden state transition mode

4.3 損傷定位和定量

仍以工況11～工況14 為例，按照式(16)計算得到各工況的損傷定位向量如表2 所示。dlv3值最大表明裂縫位于傳感器3，與實際裂縫位置相符。

表2 損傷定位結果Table 2 Damage location results

對采集的結構加速度響應數據進行二次積分和高通濾波得到結構位移響應數據，然后利用式(17)估計裂縫寬度。如表3 所示，損傷工況11～工況13 裂縫寬度的估計誤差均控制在5%左右；工況14 裂縫寬度的估計偏差較大，其原因可能是加速度積分求得的位移與結構真實位移存在的誤差與裂縫寬度相當。

表3 裂縫寬度估計結果Table 3 Crack width estimation results

5 結論

本文假定損傷會導致結構響應出現非線性特征，提出基于MSVAR 模型的非線性損傷識別方法。通過數值算例和模型試驗得出以下結論：

(1)利用MSVAR 模型的隱狀態平滑概率構造信息熵能夠準確預警非線性損傷；

(2) MSVAR 模型第一階自回歸系數矩陣的對角線元素包含裂縫位置信息；

(3)確定裂縫位置后，采用平滑概率從0 突變到100%時刻的層間位移可近似裂縫寬度；

(4)本文所提方法對寬度較小的裂縫較為敏感，但隱狀態個數取值過小會降低損傷預警精度。

需要指出的是，由于時域的時間序列模型通常易受噪聲影響，將基于MSVAR 模型的損傷識別方法應用到實際工程需要解決環境噪聲的干擾問題。另外，本文只針對裂縫類型損傷進行研究，后續需要驗證所提方法在其他類型損傷(比如鋼筋塑性變形等)的適用性。