智慧型無線感應器系統平臺建構與其在土木結構診斷應用之研究---子計畫：利用無線監測系統於土木工程之微動及破壞試驗量測與分析(I)

行政院國家科學委員會專題研究計畫 成果報告

智慧型無線感應器系統平臺建構與其在土木結構診斷應用

之研究--子計畫:利用無線監測系統於土木工程之微動及破

壞試驗量測與分析(I)

研究成果報告(完整版)

計 畫 類 別 ： 整合型 計 畫 編 號 ： NSC 97-2625-M-009-011- 執 行 期 間 ： 97 年 08 月 01 日至 98 年 10 月 31 日 執 行 單 位 ： 國立交通大學土木工程學系（所） 計 畫 主 持 人 ： 黃炯憲 計畫參與人員： 碩士班研究生-兼任助理人員：王政淵 博士班研究生-兼任助理人員：蘇威智 博士班研究生-兼任助理人員：劉靖俞 處 理 方 式 ： 本計畫可公開查詢

中 華 民 國 98 年 12 月 10 日

行政院國家科學委員會補助專題研究計畫

■ 成 果 報 告 □期中進度報告

利用無線監測系統於土木工程之微動及破壞試驗量測與分析(一)

計畫類別：□ 個別型計畫

■

整合型計畫

計畫編號： NSC97-2625-M-009-011

執行期間：97 年 8 月 1 日至 98 年 10 月 31 日

計畫主持人：黃炯憲

計畫參與人員：蘇威智、劉靖俞、王政淵

成果報告類型(依經費核定清單規定繳交)：□精簡報告

■

完整報告

本成果報告包括以下應繳交之附件：

□赴國外出差或研習心得報告一份

□赴大陸地區出差或研習心得報告一份

□出席國際學術會議心得報告及發表之論文各一份

□國際合作研究計畫國外研究報告書一份

處理方式：除產學合作研究計畫、提升產業技術及人才培育研究計畫、列管

計畫及下列情形者外，得立即公開查詢

□涉及專利或其他智慧財產權，□一年□二年後可公開查詢

執行單位：國立交通大學

中 華 民 國 98 年 11 月 30 日

目錄

目錄... I 摘 要...II ABSTRACT ... III 第一章 前言... 1 第二章 無線監測系統驗證... 4 2.1 量測系統與方法 ...4 2.1.1 試驗設備... 4 2.1.2 試驗方法... 4 2.2 量測結果驗證 ...5 第三章 建構TVARX 及結構動態特性估算... 7 3.1TVARX 建模 ...7 3.2 瞬時模態參數估算...12 3.3 瞬時模態公式之驗證...14 第四章 方法驗證、參數探討與應用... 21 4.1 參數探討 ...22 4.2 雜訊之影響...24 4.3 與其他方法比較...26 4.4 應用於振動台試驗數據分析...30 4.4.1 待測結構物介紹... 30 4.4.2 識別結果... 31 第五章 結論 ... 33 參考文獻... 34

摘 要

本研究旨在利用時變 TVARX(time-varying autoregressive with exogenous input)模型

以準確識別出結構系統之瞬時模態參數。藉由推導時變結構系統之運動方程和TVARX 模型之對等關係，證明時變系統之瞬時模態參數可透過由位移反應架構之TVARX 時變 係數直接求得。本研究亦提出以移動最小平方差法結合多項式基底函數架構TVARX 模 型之時變係數之形狀函數。傳統基底函數展開法之函數係數為常數，本方法所用基底函 數係數亦為時間函數；故僅須使用少量之基底函數即可準確近似TVARX 之係數函數。 透過模擬之數值反應驗證各種參數對識別結果之影響，並與一些已發表之識別技巧(如 遞迴識別法配合遺忘因子、傳統基底函數展開法、權重基底函數展開法)所得之分析結 果比較，證明本研究所提識別方法較優越。最後將此方法應用於處理RC 構架振動台非 線性反應資料，識別其瞬時模態特性。識別結果與觀察到之物理現象一致。 關鍵字：時變線性系統，移動最小平方差法，TVARX，系統識別，瞬時模態參數。

Abstract

This work presents an approach that accurately identifies instantaneous modal parameters of a structure using time-varying autoregressive with exogenous input (TVARX) model. By developing the equivalent relations between the equation of motion of a time-varying structural system and the TVARX model, this work proves that instantaneous modal parameters of a time-varying system can be directly estimated from the TVARX model coefficients established from displacement responses. A moving least-squares technique incorporating polynomial basis functions is adopted to approximate the coefficient functions of the TVARX model. The coefficient functions of the TVARX model are represented by polynomials having time-dependent coefficients, instead of constant coefficients as in traditional basis function expansion approaches, so that only low orders of polynomial basis functions are needed. Numerical studies are carried out to investigate the effects of parameters in the proposed approach on accurately determining instantaneous modal parameters. Numerical analyses also demonstrate that the proposed approach is superior to some published techniques (i.e., recursive technique with a forgetting factor, traditional basis function expansion approach, and weighted basis function expansion approach) in accurately estimating instantaneous modal parameters of a structure. Finally, the proposed approach is applied to process measured data for a frame specimen subjected to a series of base excitations in shaking table tests. The specimen was damaged during testing. The identified instantaneous modal parameters are consistent with observed physical phenomena.

Key words: time varying linear system, moving least squares, TVARX , system identification, instantaneous modal parameters

第一章

前言

無線監測系統為目前監測系統之發展主流之一。無線監測系統優於傳統有線監測系 統，最明顯者為佈置省時，而且監測系統每一頻道（channel）之單價較低。另一可能潛在 之優點為可以將分析之軟體植入一監測單元（sensing unit）；但以目前之省電考量、分析速 度及記憶體之限制，此優點尚無法展現出來。例如，植入利用單一測點資料建立AR 模式 或是進行損害評估（如Lynch (2007) 所提），從損害評估之實際經驗可知此不可行，除非 感應子（sensor）被佈滿整個結構系統，如同神經佈滿人全身。 本計畫為兩年計畫，其主要目的有二：（1）應用或測試本整合案所發展之無線監測系 統；（2）利用時變系統建立損害指標。第一年將利用所發展之無線監測系統於橋樑之微動 量測；並與本系所擁有之傳統微動量測系統（Tokyo Sokushin, spc51）比較所量測速度之均 方根值（rms，以三分之一八分倍頻表示）；並進一步利用小波轉換技巧識別所量測橋樑之 模態特性。另外，亦將發展一程序利用移動最小平方差法配合多項式基底建立 TVARX

（Time Varying AutoRegressive with eXogenous inputs）模式。第二年則將應用所發展無線

監測系統於振動台試驗；借用國家地震工程研究中心之大型振動台進行RC 梁柱接頭及鋼 構梁柱接頭破壞試驗。由於實驗將進行至梁柱接頭破壞；若使用傳統有線監測系統，則其 纜線可能防礙實驗之進行。然後，再利用第一年所發展之識別技巧建立實驗數據之TVARX 模式，進而發展損害指標。 有關無線監測系統於土木工程應用之發展及相關挑戰議題，可參考Lynch (2007) 所 發表之回顧文章，於此不再贅述。在本計畫本擬應用無線監測系統於橋樑微動量測，探討 無線監測系統同步量測問題與遮蔽效應。但由於研究設備費完全被刪除，本研究僅向其他 子計畫借用無線監測系統（配合加速度計）)並同時利用本系所有之有線微動量測系統(配 合速度計)，於建築物進行微動量測，驗證該無線監測系統。 為確認本整合案發展之無線監測系統之可靠性，所量測微動反應與傳統有線量測系統

(stationary process)之假設，比較量測值速度頻譜(autospectra)。若兩量測系統所得 之頻譜類似，則進一步依一般正常程序(Golden, 1991）計算速度均方根值；並且由不同量 測系統所得反應，經系統識別分析，在相同之分析程序中比對識別模態參數之差異。 由於本計畫第一年亦將發展一套建立 TVARX 模式之方法，並依所建立之 TVARX 模 式估算結構系統之瞬時模態參數。隨時間變化之瞬時模態參數應可提供結構系統是否有損 壞之訊息（此將為第二年之工作）。因此，有必要針對TVARX 建模之相關文獻進行回顧。 時變系統常常在控制、訊號處理、經濟學、腦電圖分析、地震學等領域中被發現，可

用 來 描 述 與 追 蹤 非 穩 態 過 程 之 各 種 特 性 。TVARX 或 TVARMAX （ Time Varying

AutoRegressive Moving Average with eXogenous inputs.）模式常被用於描述時變系統之輸入 與輸出關係。考慮其模型參數為時間相依。建立TVARX模式之方法可概略分為兩大類： （1）遞迴識別技巧：包括遞迴最小平方演算法（RLS）(Liung, 1987; Johansson, 1993; Parkum

et al., 1992)、遞推輔助變量演算法（RIV）( Liung, 1987)、擬線性遞迴演算法（PLR）(Liung,

1987)等方法。此些方法可以分析緩慢變化之穩態過程，但是卻無法處理急劇變化之系統。 結合卡氏濾波演算法（Kalman filtering algorithm）或貝式分析（Bayesian inference）可以 識別急劇變化過程，但卻容易過量估算模型參數以及估算之變異量過高（Ljung and Ljung, 1985； Ljung 1987；Kalouptsidis and Theodoridis,1993）。

（2）基底函數展開法：利用一系列之基底函數對每一個時變系數進行展開，然後再利用

最小平方差演算法（LS）估算此些基底函數對應之常係數（Grenier 1983, Ljung and

Soderstrom 1983, Charbonnier et al. 1987, Tsatsanis and Giannakis 1993, Rajan and Rayner

1996）。有幾種基底函數常被使用，包括傅利葉序列(Marmarelis, 1987)、Legendre Polynomials

(Zou et al., 2002)、Walsh function (Zou et al., 2002)、以及小波(Tsatsanis and Giannakis, 1993)。於吾人94年度國科會計畫中曾提出小波基底配合權重最小平方差演算法(weighted least squares approach)。此些方法所得識別結果的準確度與其所選擇之基底函數有很大的 關聯。雖然有各式各樣之基底函數可使用，但是目前並沒有可靠的方法指引選擇最佳之基 底函數及項數之多寡。使用過多項式基底函數亦通常會造成數值困難，以致無法求得準確

本研究將提出新的方法由結構系統之動態反應架構適當之TVARX模式；並依TVARX 模式中之時變系數估算該結構系統之瞬時模態參數。前述之參考文獻均侷限於TVARX模式 之建模；純粹建立TVARX模式對土木工程之應用不是那麼明顯。例如，主動控制之結構系 統於不同之輸入常導致不同之控制輸出。又例如結構非線性反應亦可用TVARX模式模擬其 該次非線性反應之輸入與輸出關係；但該輸入與輸出關係通常不適用於下次之非線性反 應。但利用TVARX模式估算結構系統之瞬時模態參數，可了解該結構系統於動態反應中之 行為；可能進一步應用於損害評估。 本研究所提之TVARX模式建模方法乃利用多項式基底展開TVARX模式中之係數函 數，再透過移動最小平方差法（moving least square）建立該係數函數之形狀函數。此方法 常用於無網格有限元素法（meshless finite element method）（Liu, 2003）。如此，各多項式 基底函數所對應之系數不是常數，而是時間函數；故本方法只須較少之基底函數即可準確 近似TVARX模型中之係數函數。

為驗證本研究所提程序之可行性及準確性，將分析一單自由度時變系統之數值模擬地 震反應，與前人所提方法之結果比較；並探討本程序中重要參數之影響，以利將來使用者 方便應用。最後亦將此識別方法應用於RC構架之振動台試驗數據分析。

第二章

無線監測系統驗證

2.1 量測系統與方法

本次實驗為使用無線監測系統(配合加速度計)，並同時利用有線微動量測系統(配合速 度計)，於交通大學工程二館進行微動量測試驗，以驗證該無線監測系統。施測方向為短向 (X 向)及長向(Y 向)，工程二館之簡易示意圖及測點概略位置如圖(2.1)所示。

2.1.1 試驗設備

試驗所使用之無線監測系統的集錄系統和感應模組如圖2.2 與圖 2.3 所示，其規格如

表2.1 所列。感應模組 CC2430/31 Sensor Board 裝配單軸加速度計、16 位元 A/D 轉換及 SPI

操作介面，其所需之供電電壓為 3.3V，最大測量範圍為±1g。 所使用有線微動量測系統(SPC-51)及速度計如圖 2.4 與圖 2.5 所示，其規格如表 2.2 及 表2.3 所示。整套系統之解析度可達_5*106_{(cm/sec)，雜訊(spectral noise)於 1~100Hz 間可} 小於2*106(cm/sec)(均方根值，rms)。值得注意的是，SPC-51 集錄系統之低通濾波器自動 設定於 3 1 取樣頻率處。故吾人只會量到頻率低於 3 200 Hz 之振動量。

2.1.2 試驗方法

本試驗於工程二館進行，選定一測點，於4 個樓層中分別擺設兩個 VSE15D 速度計

以及一個CC2430/31 Sensor Board 感應模組。各個樓層之 VSE15D 速度計分佈於短向(X 向)

及長向(Y 向)以量測相對應之反應。由於祇是四個 CC2430/31 Sensor Board 感應模組，是 故短向(X 向)與長向(Y 向)分兩次進行量測。兩種感應計擺設方位如圖 2.6，現地量測照片

施測方式說明如下： 1.量測時間： 每段量測時間皆紀錄5 分鐘，取樣頻率為 100Hz，故每個測點共取 30000 點。 2.量測振幅範圍 兩種感應計分別量測結構之速度反應及加速度反應。因微動試驗之環境外力振幅極 小，故有線微動量測系統所設之振幅範圍為±1kine，無線監測系統則設為±1g。

2.2 量測結果驗證

為驗證無線監測系統之可靠性，所量測微動反應與傳統有線量測系統所得者比較。由 於一為量測加速度反應，另一為量測速度反應；故先利用平穩過程(stationary process)之假 設，以比較量測值速度頻譜。圖2.8 和圖 2.9 為工程二館各樓層之無線及有線量測於 X 向(短 向)和 Y 向(長向)的微動量測試驗歷時，由圖 2.8 及圖 2.9 可初步看出，無線監測系統之量 測結果並不理想。 進一步分析兩種量測數據，透過下式估算各測點之速度及加速度單邊能譜密度函數

 

_





  nc k k d x x V f T T n f G 1 2 , 2   (2.1)

 





2 1 2 2 , c n k k d xx A f T n T G f         



 (2.2) 其中：V_k



f,T



為第k 段速度反應歷時之 Fourier transform、A f T_k



,



為第 k 段加速度反應

線微動量測系統與無線監測系統所得之結果，由無線監測系統所得加速度頻譜透過

 

f f G

 

f G_xx ) 2 _{xx} 2 (    , 轉換成速度頻譜。由有線微動量測系統與無線監測系統所得之速度頻譜如圖2.10 及圖 2.11 所示。 由圖 2.10 及圖 2.11 可發現，有線微動量測系統與無線監測系統所得頻譜非常的不一 樣，值之大小與圖形均有非常大之差異。有線微動量測系統於短向(X 向)在不同測點之分 析結果，可清楚看出約在5Hz 附近有明顯之峰值；長向(Y 向)則是在 4Hz 附近有明顯之峰 值。而無線監測系統之分析結果，雖然均有明顯峰值，但仔細觀察可發現，其峰值成規律 性出現，約每5Hz 出現一峰值，造成此現象有可能是 CC2430/31 Sensor Board 電路不穩所 造成之量測誤差或者是轉檔之編譯程式撰寫不佳所致。

第三章 建構

TVARX 及結構動態特性估算

3.1 TVARX 建模

在一時變線性系統之反應，時變ARX（TV-ARX）模式可表示為

 

t

   

t t i J

  

t t j



 

t j j I i i y Θ f an Φ y 



 



   1 0 . (3.1) 其中，y(ti)是在tit時刻下量得之反應向量；f(ti)是在tit量得之外力輸入向量； ) (t i Φ 和Θ_j(t)為待測之時變係數矩陣；I與J則分別代表所考慮模型之階數；a_n(t)為誤差 向量，由於量測噪訊或模型誤差所造成。TVARX 通常被應用於架構時變線性系統或非線 性系統之輸入與輸出關係。當利用TVARX 模式建構一非線性系統之輸入及輸出時；即是 將與位移或速度反應相關之系統特性(勁度或阻尼，為時間之隱函數)以一等價之時間顯函 數替代。欲從量測動態反應建立TVARX 模式，吾人將提出一利用多項式基底函數配合移 動最小平方差法決定式(3.1)之時變係數矩陣。 TVARX 模型中之係數函數以多項式基底函數展開。將式(3.1)之時變係數矩陣利用多 項式基底函數展開得

 

i i Ni n n n i i t A t A P Φ ~ 0 ˆ ˆ ˆ  



 ， (3.2a)

 

j j N n n n j j j t t B B T Θ ~ 0 ˆ ˆ ˆ  



 ， (3.2b) 其中，

     

 





T I I I I P Ni i t t  t 2  ， (3.3a)

     

 





T I I I I T Nj j  t t2  t ， (3.3b)



i i iNi



i A A A A~  _{0 1 } ， (3.3c)



j j jNj



j B B B B~  _{0 1 } ， (3.3d) n i ˆ A 與B 為待定係數矩陣，_{j ˆn}

 

I 為維度與係數矩陣一致之單位矩陣。

式(3.2)中之待定係數A 與_iˆn B 以權重最小平方差法（_{j ˆn} weighted least square）進行求解。

令Φ_i

 

t_k 與Θ_j

 

t_k 之真值分別為Φik與Θ 。待定係數可透過式(3.2)所得之jk Φi

 

tk 及Θj

 

tk 分別與其真值之誤差函數最小化求得。以式(3.2a)為例，將求A~_i之誤差函數定義為

 

_













    li l il i i l n il i i n t t t 1 ~ , ~ Φ P A W Φ P A E T ， (3.4) 其中：









 

I W t_n~,t_l wt~_n,t_l ； l_i為描述Φ_i

 

t 上所取之節點總數，其值遠小於反應總取樣點數；



tn tl



w , 為權重函數，此權重函數亦可視為一窗函數，須為一正定函數(positive definite function)。本研究中選用指數型權重函數：





          1 0 1 , 2 ~ ~ d d e t t w _{n m} d ， (3.5) 其中，d  t_n~t_m~ d_m；d 代表_m w支撐長度。此函數圖形如圖3.1 所示；其中 0.3。 對函數E 最小化可得

 

0 AE    i ~ ， (3.6) 將誤差函數代入上式；整理可得

 

t _{i i}

 

t i iA ΦQ A~  ， (3.7) 其中

 

_

  

  

  li l l i l n l i i t t t t t 1 , P W P A T ， Q

 



Q Q Q



T i l i i i i t  1 2  ， Qil Pi

  

tl W tn,tl



，



i i _il



i Φ Φ Φ Φ  _{1 2 } ， (3.8) 則式中A~_i之解為

   

t _i t i i i 1 ~ _{ΦQ A} A ， (3.9) 上式所求之A~_i為時間t 之函數；將之代入式(3.7)得

 

t _{i pi}

 

t i , ~ Φ Φ Φ  ， (3.10)

 

t i p, ~ Φ 為有限元素中所稱之形狀函數（shape function）矩陣

 

     

T P A Q Φ_pi t _i t _i1 t _i t , ~ _  _{， (3.11)} 同理，Θ_j

 

t 可用相同方式展開

 

t _{j p j}

 

t j , ~ Θ Θ Θ  ， (3.12)



j j jlj



j Θ 1 Θ 2 Θ ˆ Θ   ，Θp j

 

t Uj

     

t Bj t Tj t T 1 , ~ _  _，

 

  

  

l j l n l l l j j t t t t t j T W T B T _, ˆ 1



  ，

 





T U U U U j l j j j j t  1 2  ˆ ，Ujl Tj

  

tl W tn,tl



， (3.13) 而lˆj為描述Θj

 

t 上所取之節點總數。如果對TVARX 模型中之各個時變係數取相同之結點 數與相同之多項式基底，則簡化以上表示式。 在式（3.10）與式（3.12）中未知數為Φ 和_i Θ_j，以最小平方差法求解。其誤差函數之 定義為

 







 





  N n n n T n n t t 1 a a E ， (3.14) 其中，N 為資料點數。由式(3.1)可知誤差可由下式表示

   

   

  



 

   

  



 

_                              













      J j j t j I i i ti J j j p j I i i pi J j j I i i n t j t t i t t t j t t i t t t t 0 , 1 , 0 , 1 , 0 1 ~ ~ Ω Θ Γ Φ y f Θ Θ y Φ Φ y f Θ y Φ y a ， (3.15) 其中

 

 

 



t i t i t i



i l i p i p p,i, i t  Φ y  Φ y  Φ y  Γ, 1 ,,2 _,, ~ ~ ~  ， (3.16a)















t j t j t j



j l j p j p j p j t  Θ f  Θ f  Θ f  Ω , , ,1 , ,2 _,,ˆ ~ ~ ~  ， (316b)

式(3.14)之誤差最小，可得



  

y Φ Γ Θ Ω

 

0 y Φ Φ E _                    







   N n J j j t j I i i ti n n i i t i t 1 1 , 0 , ~ 2 ， (3.17a)



  

y ΦΓ Θ Ω

 

0 f Θ Θ E                     







   N n J j j t j I i i ti n n j j t j t 1 1 , 0 , ~ 2 ， (3.17b) 將以上之聯立方程組以矩陣形式表式可得 m m m mF V V C VT  T ， (3.18) 其中                J t t t I t t t J t t t I t t t J t t t I t t t m n n n n n n,1 ,2 , ,0 ,1 , , 1 , 0 , , 2 , 1 , , 1 , 0 , , 2 , 1 , 2 2 2 2 2 2 1 1 1 1 1 1 Ω Ω Ω Γ Γ Γ Ω Ω Ω Γ Γ Γ Ω Ω Ω Γ Γ Γ V               ， (3,19a)

   

 





T n t t t y y y F _{1 2 } ， (3.19b)



_T



T J T T T I T T m Φ Φ Φ Θ Θ Θ C  _{1 2  0 1 } ， (3.19c) 對 (3.18)式求解即可得到每個節點之近似真值，由此即可得到此系統於某個時刻t  下之t_n 時變係數Φ_i

  

t_n i1,2_I



。 架構 TVARX 之各個時變係數時，可考慮不同之多項式階數N 與_i N_j以及不同之節點 數 l 與_i lˆ 。 於 實 際 應 用 時 若 取 相 同 之 多 項 式 階 數 以 及 相 同 之 節 點 數 ， 令_j          ˆ ˆ   ˆ 

係數可使用相同之形狀函數，將可使運算過程更有效率。

3.2 瞬時模態參數估算

從量測反應建構適當之TVARX 模式，通常吾人欲瞭解該時變線性結構系統之動態特 性隨時間變化之行為。此等信息有助於判斷結構系統之損傷情形。在瞬時t 下，Φ_i

 

t 與

 

t j Θ 均為常數矩陣，故在瞬時t 下，TVARX 模式即對等於非時變之 ARX 模式。依非時 變ARX 模式估算動態特性之方法，令





 

 

 

 

 

                   t t t t t t I I I I 1 2 3 1 0 0 0 0 0 0 0 0 0 0 0 0 ) ( Φ Φ Φ Φ Φ I I I G           (3.20) 在時變系統

 

G 將隨時間而改變。

 

G 之特徵值及向量是與結構系統之瞬時動態特性有直接 關係。如同非時變系統之推導可知

 

G Ψk kΨk ~ ~ _ (3.21) 其中Ψ~_k及_k為

 

G 之第k 特徵向量及特徵值。在時變系統中，Ψ~_k及_k亦是時間之函數。 令                ) ( ) 2 ( ) 1 ( ~ ~ ~ ~ I k k k k Ψ Ψ Ψ Ψ 

其中，Ψ~(_ki)為一

 

l1 之向量，l為量測自由度。由於

 

G 之特殊構造，可得 ) 1 ( ) ( ~ ~ _ i k k i k Ψ Ψ  ( i= 1,2,…,I ) (3.22a) 且 ) ( ) 1 ( 1 0 ~ ~ ~ I k k j k I j j n Ψ Ψ φ     



(3.22b) 從上式兩關係式（（3.22a）及（3.22b））可得                  1 (1) ) 1 ( 2 ) 1 ( ) 1 ( ~ ~ ~ ~ ~ k I k k k k k k k λ λ λ Ψ Ψ Ψ Ψ Ψ  (3.23) 及Ψ 代表量測自由度之第(_k1) k 模態。 式（3.20）中之特徵值常為複數，成雙或對之共軛根。令_k a_k ib_k，則結構系統 之瞬時擬自然振動頻率及阻尼比為 2 2 ~ k k k      ， (3.24) 其中： k k k a b t 1 tan 1     ， (3.25) t  為時間增量(即取樣頻率之倒數)。

3.3 瞬時模態公式之驗證

依上節所述程序所得瞬時擬自然振動頻率及阻尼比之正確性，本節以一單自由度系統 透過嚴謹數學證明之。單自由度時變線性系統之運動方程式為：

           

t x t c t x t k t x t f

 

t m      (3.26) 其中，m，c，k分別為該系統之質量，阻尼係數及勁度，為時間之函數；f

 

t 為外力；x

 

t ，

 

t x 及x

 

t 分別為位移，速度及加速度反應。一時變系統之瞬時模態參數可透過式(3.26) 表示，其定義方法類似於線性非時變系統之模態參數 ) ( ) ( ) ( 2 ) ( t m t k t f t _n n     ； ) ( ) ( 2 ) ( ) ( t t m t c t n    ， (3.27) 其中，_n(t)與(t)分別為瞬時自然振動頻率與阻尼比，為時間函數。 單自由度系統TVARX 模式為

 

t

   

t yt i

  

t f t j



a

 

t y J _n j j I i i     





 1 0   ， (3.28) 在以下討論TVARX 與運動方程式對等關係時，不考慮 TVARX 中之量測誤差項。

吾人可利用中央差分法（central difference approach）將式(3.26)離散化得

 

 

   

 

 

 

 

 

 

 

   

t x t t t c t t m t x t t m t k t f t t x t t c t t m                                  2 2 2 2 2 2 ， (3.29) 其中，t為時間增量。為對等式(3.29)與式(3.26)將式(3.29)表示成：

 

t  ₁

   

t xt1  ₂

  

t x t2



 ₁

   

t f t1 x    ， (3.30) 其中，x

 

ti 即代表x



tit



，

 







  



t t



t t t m t t k t            ₁ 2 2 ，

 



  

_



_

  

t t t t t c t t t m t             ₂ 2 2 ，

 

_

_

t t t     ₁ 1 ，









 



 

t



t t c t t t m t t           2 2  。 由於現地量測時，通常只量測速度或加速度反應，故若利用速度或加速度直接建立 TVARX 模式。將式(3.26)分別對 t 微分一次及二次得：

   

t v t



m

   

t ct



v

     

t



ct k t



v

     

t k t x t f

 

t m _  _  _  _      ， (3.31a)

   

t a t



m

   

t c t



a

 

t



m

 

t c

   

t k t



a

   

t



c t k

 

t



v

     

t k t xt f

 

t m _  2_  _  _ 2_   _ 2     ， (3.31b) 其中，v

 

t 與a

 

t 分別代表速度與加速度反應。當欲用速度反應架構TVARX 時，由式(3.26) 得

             



f t mt vt c t v t



t k t x  1    ， (3.32) 代入式（3.31a）得

           

t v t c t vt k t v t f

 

t m      ， (3.33) 其中

 

   

   

_{ }

t k t k t m t c t m t c      ，

         

_{ }

t k t k t c t c t k t k     ，

 

   

_{   }

f t t k t k t f t f     ， 式(3.33)可應用中央差分法表示為

 

t  ₁

   

t v t1  ₂

  

t v t2



 ₀

   

t f t  ₁

   

t f t1  ₂

  

t f t2



v      ， (3.34) 其中

 





_



_

  

t t t t t m t t k t            ₁ 2 2 ，

 



  

_



_

  

t t t t t c t t t m t             ₂ 2 2 ，

     

t  t  tt  2 1 0 ，

 







t t

 

t t



k t t k t          ₁  ，

     

t t tt     2 1 2 ，









 



 

t



t t c t t t m t t           2 2  ， 當欲用加速度反應架構TVARX 時，由式(3.33)得

 

_{           }



f t mt at c t at



t k t v  1    ， (3.35) 將式(3.32)與(3.35)代入式(3.34)，整理得

           

t at c t at k t a t f

 

t m _ ~ _ ~  ~ ， (3.36) 其中

 

   

   

_{ }

t k t k t m t c t m t c  _    ~ _，

         

_{ }

t k t k t c t c t k t k~     ，

 

   

_{   }

f t t k t k t f t f     ~ ， 同樣地，應用中央差分法至式(3.36)可得

 

t  ~₁

   

t a t1  ~₂

  

t a t2



 ~₀

   

t f t  ~₁

   

t f t1  ~₂

  

t f t2



a      ， (3.37) 其中，

 





_



_

  

t t t t t m t t k t            ~₁ ~ _~2 2 ，

 



  

_



_

  

t t t t t c t t t m t             ~_{2 ~}2 ~ 2 ，

 

_

_{  }

_{  }





_

_{  }





_

_                       t t k t t t k t t k t t t k t t t t 2 2 1 ~ 1 ~ 2 0     ，

 

_

_{  }



 























 

 





                               t t k t t k t t k t t k t t k t t k t t k t t k t t t t  ₂   2 2 1 2 ~ 1   ，

 

_

_{  }

_{  }





_

_{  }





_

_                       t t k t t t k t t k t t t k t t t t 2 2 1 ~ 1 ~ 2 2     ，









 



 

t



t t c t t t m t t           2 ~ ~ 2  ， 因此，單自由度之線性時變結構系統，若以位移反應來建構TVARX 模式；則其理論 階數

   

I,J  2,1 ，且₀

 

t 0。若以速度或加速度反應來建構 TVARX 模式，則其理論階 數

   

I,J  2,2 。 當用位移反應建立TVARX 模式時，得如式(3.30)所示 TVARX 模式，代入式(3.20)得

 

G _0 1_ (3.38)

令其特徵值為_1,2 a₁ib₁。其中 t c m k t m a      2 ) ( 2 2 1 ；





t c m c k t mk t b        2 ) ( 4 2 2 2 1/2 1 (3.39) 上式及以下各方程式之推導，為了簡化該表示式，省略m, c, k 函數中之(tt)。代入 式 (3.25)，可得 ) 2 2 ln( 2 1 1 t c m t c m t        ， (3.40a)





_] ) ( 2 ) ( 4 [ tan 1 2 2 / 1 2 2 2 1 1 t k m c t k mk t t           ， (3.40b)

利用Taylor’s expansion，將上式之 ln 項及_tan1_{項展開可得}

                  2 1 2 T t O T   ， (3.41a)

 

 

 

 

 

                                                    4 2 2 3 2 3 2 1 2 3 2 1 2 2 2 1 2 1 1 16 1 24 1 2 3 ! 3 1 1 2 T t O T t T T T          ， (3.41b) 其中 k m T 2 ， m c T   4  ，







n



T t O  代表含t T階數大於（或等於）n 之項。 因為在利用有限差分法離散化運動方程式時，要求 Tt 1；因此，上式之₁及₁可簡 化成 2 / 1 2 1 1 (1 ) 2 , 2 _  _     T T ， (3.42) 故，依式(3.24)所得瞬時擬自然振動頻率及阻尼比為 m k T    ₁ 2 ；₁ ， (3.43) 因此，依前節所得之擬瞬時模態估算結果即為系統在 tt時刻之瞬時自然振動頻率及阻 尼比。 若是利用速度或加速度建立TVARX 模式，再依上節所述估算瞬時擬自然振動頻率及 阻尼比，則分別得速度反應建模之估算結果 m k v  1  ； m c T v _  4 1  ， (3.44) 其中：

   

   

_{ }

t k t k t m t c t c    ，

         

_{ }

t k t k t c t c t k t k   _{ }  ， 以及加速度反應建模之估算結果 m k a ~ 1   ； m c T a ~ 4 1 _   ， (3.45) 其中

   

   

_{ }

t k t k t m t c t c    ~ _，

_{         }

 

t k t k t c t c t k t k~      ， 由以上關係式可明顯看出，由位移、速度與加速度反應估算結構系統之瞬時模態 反應，此三組反應所得之估算結果並不一致。由位移反應所得之瞬時模態與理論值完全一 致；而由速度反應估算之結果將包含c

 

t 與k 所構成之誤差項；由加速度反應估算之結

 

t 果將包含c

 

t 、c

 

t 、k 與

 

t k 所構成之誤差項。

 

t

第四章

方法驗證、參數探討與應用

為驗證此識別方法對時變系統之識別效果，將透過一系列之單自由度數值模擬反應對 此識別流程進行驗證，並且進行各種參數之探討以掌握此識別方法之特性。另外，並與使 用遞迴最小平方差法(recursive least squares approach, Ljung, 1987))及加權基底函數法 （Weighted Basis Function approach, WBF, Niedzwiecki, 2000）

考慮一單自由度時變運動方程式：

 

t 2

     

t t x t

   

t 2x t a (t) x   _n  _n  _g (4.1) 其中探討四種不同型式之瞬時模態參數定義如下： Case 1：平緩變化系統：               t t n 60 2 sin 5 . 0 5 . 1 2 ) (    ;         t t 60 2 sin 2 4 ) (   ， (4.2) Case 2：週期變化系統：               t t n ₁₀ 2 sin 5 . 0 0 . 1 2 ) (    ;         t t 10 2 sin 5 . 2 5 ) (   ， (4.3) Case 3：跳躍變化系統：                 30 15 5 . 7 15 0 5 ) ( 30 15 5 . 0 15 0 1 ) ( t t t t t t n   ， (4.4) Case 4：折線變化系統：

























                                                                     30 27 5 27 23 5 . 2 19 4 5 . 2 2.5 19 23 19 11 5 . 7 11 4 5 . 2 7.5 7 11 7 3 5 3 4 5 . 2 5 0 3 ) ( 30 27 27 23 23 19 19 11 11 7 7 3 3 0 1 5 . 1 19 4 5 . 0 1.5 5 . 0 11 4 5 . 0 0.5 1 3 4 5 . 0 1 ) ( t t t t t t t t t t t t t t t t t t t t t t n   ， (4.5) 相較於平緩變化系統，週期變化系統之瞬時模態參數變化較為劇烈；而跳躍變化系統 則於某一時刻下為不連續之變化，其瞬時模態變化劇烈程度更甚於週期變化系統；雖然折 線變化系統之瞬時模態參數為連續變化，但是其變化曲線為一次微分不連續。 以 Runge-Kutta 求出此時變系統受地震輸入(a_g(t))之動態反應，取時間增量(t)為 0.001 秒。輸入之地震反應及其頻譜圖如圖 4.1 所示。以下將利用位移反應進行識別。而平 緩變化、週期變化、跳躍變化以及折線變化系統之時間歷時反應及其頻譜分別示於圖4.2

4.1

參數探討

利用位移反應建構TVARX 模型，將 TVARX 之時變係數以移動最小平方差法建構各 個節點之形狀函數展開，再求得各形狀函數所對應之節點值。分析中主要探討此方法之各 項參數（如多項式基底階數N 、加權函數之支撐長度_p d 與節點數_m l ）對識別結果之影響。_n 於TVARX 模式與運動方程之對應關係之推導中得知，TVARX 之理論階數為

   

I,J  2,1 ； 故對無雜訊干擾系統進行分析時模型階數將取

   

I,J  2,1。以下將探討多項式階數對各組 時變系統識別之影響，而識別結果與實際時變特性之相對誤差定義如下：

100% ) ( ) ( ) (   i true i true i id t t t    ， (4.6) 其中_true與_id分別代表瞬時自然振動頻率或阻尼比之實際值與識別所得之結果。而相對誤 差之平均值與標準差分別記為與 。 以多項式基底函數對時變係數展開時，選用固定之多項式階數N_p 2，圖 4.3 為各組 時變系統比較不同支撐參數在考慮不同結點數下識別所得之相對誤差平均值以及標準 差。分析中選用之加權函數支撐長度分別為d_m 2、d_m 4以及d_m 6。由圖4.2 可發現： (a)識別平緩變化系統之瞬時自然振動頻率時，當d_m 2時，取節點數l_n 35將有較好 之結果，當d_m 4與d_m 6時，於圖中選用之節點數範圍內其識別之平均誤差均遠小於 1%；在瞬時阻尼比之識別時，當d_m 2時，取節點數l_n 35將有較好之結果，當d_m 4與 6  m d 時，於圖中選用之結點數範圍內其識別之平均誤差均遠小於20%。 (b)識別週期變化系統之瞬時模態時，其識別之平均誤差之趨勢與平緩變化系統相似， 唯有在d_m 4與d_m 6時當節點數l_n 30時，其識別之平均誤差略高。 (c)識別對跳躍變化系統之瞬時模態時，當d_m 2時，節點數l_n 35其頻率誤差之平均 介於1%~2%之間，而阻尼比誤差之平均介於 10%~20%之間；當d_m 4與d_m 6時，其頻 率與阻尼比之識別誤差分別無法穩定控制於2%與 20%以下。 (d)識別折線變化系統之瞬時模態時，其識別結果之平均誤差趨勢與平緩變化及週期變 化系統相似。 當取d_m 2時，圖 4.4 為各組時變系統比較不同多項式階數在考慮不同結點數下識別 所得之相對誤差平均值以及標準差；分析中選用之多項式階數分別為N_p 1、N_p 2以及

(a)識別平緩變化系統之瞬時模態時，當N_p 1以及N_p 2時，於節點數l_n 35有較 小之識別誤差；當N_p 3時，於選用之節點數範圍內其頻率與阻尼比識別之平均誤差分別 遠小於1%與 20%。 (b)識別週期變化系統之瞬時模態時，其識別誤差之趨勢與平緩變化系統相似，唯在 3  p N 且當節點數為l_n 20時，其頻率識別之平均誤差略高於1%。 (c)識別對跳躍變化系統之瞬時模態時，此三組多項式階數於節點l_n 35時，其識別頻 率之平均誤差介於1%~2%之間，而識別阻尼比之平均誤差介於 10%~20%之間。 (d)識別折線變化系統之瞬時模態時，其識別之平均誤差趨勢與平緩變化及週期變化系 統相似。圖4.5 所示為取N_p 2以及l_n 40所得之識別結果。 以傳統多項式基底函數展開法架構TVARX 模式並識別系統之瞬時模態特性，須引入 較高階基底函數方可能得到較好之識別結果；但使用高階之多項式基底時，易造成病態矩 陣並導致數值困難。以加權基底函數法建模並進行瞬時模態之識別，雖然可使用較少之基 底函數以獲得不錯之結果，但須於每一時刻計算基底展開之係數，分析過程相當耗時。以 移動最小平方差法建模並識別，不僅可以少量之基底函數對時變係數展開，且計算形狀函 數對應之係數僅須進行一次反算，其識別效率優於傳統基底函數展開法與加權基底函數 法。最後，值得一提的是，利用多項式基底函數描述不連續函數（如case 3），其效果如所 預期地不理想。

4.2 雜訊之影響

為瞭解雜訊對識別品質之影響，於數值分析所得case 2 之相對位移反應及輸入反應中

加入白色雜訊，使得NSR=5％。然後，對此含有雜訊之訊號，進行識別，並與未含雜訊之

分析結果比較。以下將選用固定之多項式N_p 2與節點數l_n 40，探討當權重函數之支撐

參數分別為d_m 2、d_m 4及d_m 6時，對應不同TVARX 模型階數

 

I,J 之識別結果。

在求得模型中之係數後，欲進一步評估此模型是否能準確描述結構系統之輸入與輸出

關係，以下引入AIC（Akaike’s Information Criterion, Akaike, 1973）及 FPE（Final Predication

Error, Akaike, 1970）等兩種評估準則。通常利用前一步預測值（one-step ahead prediction） 與實際量測值之誤差定義損失函數，即

 



  N k k t V 1 2 n a (4.7) 一般而言，模型結構越複雜，其損失函數之值當然越小。因此，必須導入關於模型之複雜 度（即組成模型之參數），此兩種評估準則分述如下： 1. AIC 是用以量測最似然推測法所得模型之不吻合度。其定義為：



max. probability

 

2 parameter



ln 2    AIC ， 此時，AIC 值越小代表模型越接近真實系統。於計算 AIC 值時，或然率之計算是必要 的，但是在預測誤差為正規性之情況下，其表示法可改寫為：                V N n AIC _{ln 1} 2 u _{， (4.8)} 2. FPE 之定義式如下： V N N n N n FPE u u _{ }    1 / 1 1 ， (4.9)

圖 4.6 與圖 4.7 分別為利用不同d 值結合不同_m

 

I,J 所得之AIC與FPE 值變化曲線。

考慮不同d 值，均在_m

 

I,J 18時有最小AIC值。不同d 值，均在_m

 

I,J 15時均有最小

FPE 值。而此三組d 值所得之_m AIC與FPE 值變化相當一致。比較此兩種指標之分析結果

可知，以AIC值做為指標所決定之模型階數略高，但基本上兩者所決定之模型階數差異不

大。

考慮支撐參數與

 

I,J 對識別所得瞬時模態參數相對誤差之影響。藉由相對誤差之變化

曲線與AIC或FPE 值相互參照，以瞭解藉由AIC或FPE 值與識別誤差之關聯性。圖 4.8

為各組支撐參數於不同

 

I,J 模態識別誤差之平均值與標準差。由瞬時自然振動頻率識別誤 差之平均值與標準差變化曲線可看出，其誤差平均值約在

 

I,J 20~25有最小值；標準差 約在

 

I,J 20之後趨於穩定。此外，在

 

I,J 15之後識別誤差之平均值與標準差均在2% 以下。隨著

 

I,J 增加，阻尼比識別誤差越小；約在

 

I,J 17之後，阻尼比識別誤差之遞 減曲線逐漸趨於平緩。此外，支撐參數為d_m 4及d_m 6時，於

 

I,J 16之後其識別誤差 之平均值與標準差均在20%以下。值得一提的是，以支撐參數d_m 2時有最大之識別誤

差，此結果與AIC與FPE 值之變化趨勢相反。但基本上，以AIC與FPE 值所決定之最佳

模型階數與識別誤差之分析結果相當契合。 經由以上之分析，選擇以支撐參數d_m 4，模型階數

 

I,J 20進行識別。圖 4.9 為 case 2 於輸入與輸出訊號加入雜訊之識別結果。識別結果顯示，在訊號中引入雜訊仍能準 確識別瞬時自然振動頻率之變化情形；雖然對於阻尼比之識別結果較不理想，但是其誤差 尚在可接受範圍。

4.3 與其他方法比較

遞迴識別法亦常被利用於識別隨時間改變之參數。其優點為運算速度快且可線上 (on-line)執行；其缺點為追蹤參數隨時間改變之準確度較差，尤其在分析具雜訊之反應。 遞迴識別法分析步驟可歸納如下（Ljung, 1987）：

  

ˆ 1

    



  

ˆ 1





ˆ_{t  t L t y t }T _{t t }  (4.10)

 

_{ }



_{  }

  

_{  }

t t p t t t t p t L _T     1 1     (4.11)

 

_{  }





_{ }

     

_{  }

_{  }



_             t t p t t t p t t t p t p t t p _T T       1 1 1 1 1 (4.12) 上式中，y

 

t 為該時刻下之反應；

 

t 為該時刻以前之輸入與輸出組成之向量；ˆ 為代

 

t 測係數，

 

t 為遺忘因子。分析中，取

 

t 0.95。識別過程中需要給定ˆ

 

0 與p

 

0 兩項 初始條件，而初始條件之決定可先取分析資料最前段之100 點以最小平方差法進行估算。 權重基底函數展開法主要是為解決使用傳統基底函數展開法中需使用大量基底函數之問 題。其原理簡述如下（Niedźwiecki, 2000）: 對式(3.1)中之 TVARX 係數函數以多項式函數(可用其他基底函數)展開，即令：

 

_

  Ni n n n i i t t 0 ˆ ˆ ˆ A Φ ，

 



  Nj n n n j j t t 0 ˆ ˆ ˆ B Θ ， (4.13) 其中，A 和_{i ˆn} B 為待測之係數矩陣。可利用加權最小平方差法求解_{j ˆn} A 和_{i ˆn} B 。定義誤差估_{j ˆn} 計函數為

 

_



    





  M m m n T m n m n n wt t t t t ~ 0 ~ ~ ~ ~ ~ ~ , a a E ， (4.14) 其中，M~ 為估算t_n~時刻係數時，權重函數支撐長度(support)內之資料點數。 式(4.13)代入式(3.1)則可得：

 

_{ }

_{ }

                       J j j t N n n n j I i i t Ni n n n i t m m j m m t t t 1 ˆ 0 ˆ ˆ 1 ˆ 0 ˆ ˆ ~ y _~ A y _~ B f _~ a_n ， (4.15) 再代入式(4.14)；當誤差估計值最小時，誤差估計函數對各個係數取導數須為零：

 

0 AE    n i ˆ ；

 

0 BE    n j ˆ ， (4.16) 依上式整理可得：





A y B f y yT _

 

0                        



 

 

        M m i t n m t J j j t N n n m n j I i i t Ni n n m n i m n m m m j m t t t t t w ~ 0 ~ ˆ 1 ˆ 0 ˆ ~ ˆ 1 ˆ 0 ˆ ~ ˆ ~ ~, _{~ ~ ~} ， (4.17a)





A y B f y fT _

 

0                        



 

 

        M m j t n m t J j j t N n n m n j I i i t Ni n n m n i m n m m m j m t t t t t w ~ 0 ˆ 1 ˆ 0 ˆ ˆ 1 ˆ 0 ˆ ˆ , ， (4.17b) 式（4.17）可以矩陣形式簡化表示為： WF V C WV V T w w w T w  ， (4.18) 其中，



_T



T J T T T I T T w A A A B B B C  ~₁ ~_{2 } ~ ~₀ ~_{1 } ~ ，



i i iNi



i A A A A~  _{0 1 } ,





j N j j j j B B B B~  _{0 1 } ，



_{t t t}



T M~ 1 0 y y y F _ ，



w w w_M



diag _{0 1}  ~  W , w_m~ w



t_n~,t_m~



，

 







 





 





 



N t





N

 

t





N

 

t J



I t N t N t N w t t t t t t J I        f Π f Π f Π y Π y Π y Π V   1 0 2 1 1 0 2 1 ，

 









 





 





 



_T



T i t M N T i t N T i t N i t Ni t y   Π i t0 y0 Π i t1 y1 Π i t~ y M~ Π  ，

 









 





 





 



_T



T j t M N T j t N T j t N j t Nj t f   Π j t0 f0 Π j t1 f1 Π j t~ fM~ Π  ，

 



   

 

_N



T m m m Ni t~ I It~ It~i Π  _ ，

 

I 為維度與係數矩陣一致之單位矩陣。 對式(4.18)求解即可得到每一時刻下時變係數以多項式基底展開所對應之係數。由此即可 得到此TVARX 模型於某個時刻t 下之時變係數t_n~ Φ_i

  

t_n~ i1,2I



。 由式(4.14)可知『加權最小平方差法』可視為在欲分析之時間點周圍取一窗函數，且此 窗函數在靠近欲分析之時間點時有較大之權值。如此，所架構之矩陣C 為隨_w t 而改變之係 數矩陣；此外，分析過程不須太多之基底函數，即可得相當好之結果。然而，每計算一個 時間點之模態則須重新以加權最小平方差法計算係數矩陣C ，此分析過程相當耗時。當_w 所選用之權重函數為一常數（如，w



t~_n,t_m~



1），則每一時刻下，各基底函數所對應之係數 即為常數；則其所得之結果與傳統基底函數展開法一致。 圖4.10 和 4.11 比較本研究所提方法與上述遞迴識別法識別 case 1 與 case 2 之結果，其 中不考慮雜訊之影響。本研究所得結果乃利用兩階多項式基底函數配合d_m 4以及 30  n l 。該圖顯示本研究所得結果遠比遞迴識別法所得結果準確；尤其在阻尼的識別結。 在case 1 之平緩變化系統，遞迴識別法尚能提供合理之結果；但在 case 2 之週期變化系統， 遞迴識別法無法合理追蹤參數隨時間改變之特性。 利用上述傳統基底函數展開法(BF)與權重基底函數展開法(WBF)分析 case 2 之無雜訊 輸入與輸出訊號，識別之 f_n and 示於圖 4.12。傳統基底函數展開法所得之結果乃利用 20 個多項式基底函數(即{1, t, _t2 _{, ,t}19  })。權重基底函數展開法所得之結果乃利用兩階 多項式基底函數配合d_m 4。比較圖 4.12 與圖4.11 中本研究所得者可發現本研究所提方 法與權重基底函數展開法所得結果均比傳統基底函數展開法所得者精準許多。於傳統基底