自我迴歸脈衝響應函數法系統鑑定之研究
及其在修正不確定性分析模型之應用
計劃編號 : NSC 89-2611-E-002-007
執行期限 : 八十八年八月一日至八十九年七月三十一日
計劃主持人 : 洪振發 計劃參與人員 : 柯文俊、戴志豪、彭彥惇
國立台灣大學造船暨海洋工程學研究所
一、中文摘要 本計畫共分為兩大部份。第一部份為整合自我迴歸 外變數（AutoRegressive with eXogeneous variable, ARX ） 模 型 與 特 徵 系 統 識 別 運 算 (Eigensystem Realization Algorithm, ERA) 而 完 成 一 套 修 正 型 ARX/ERA 方法。此法可以由結構振動系統所量得 的輸入外力與輸出動態響應資料，建立 ARX 模型， 以求出離散馬可夫參數（Markov parameters, MP）， 再經由 ERA 法根據此參數鑑定出系統的狀態方程 式。ARX 模搭配 ERA 後，可藉由奇值分解將雜訊 去除，而將系統降階。此一優點可以藉由 ARX 模 式直接由量測資料取得馬可夫參數，藉由 ERA 免 去一般由 ARX 模型識別模態參數需要降階的困 擾。第二部份承第一部份之結果所鑑定出之模態參 數，使用特徵靈敏度反式法(Inverse Eigensensitivity Method, IEM)進一步應用至分析模型之更新修正。 本研究以一六自由度質量塊系統來說明其應用過 程。 關鍵詞：系統鑑定、自我迴歸外變數模型、特徵系 統識別運算、特徵靈敏度反式，模型更新Abstr act
This project contains two parts. In first part, we integrated the AutoRegressive with eXogeneous (ARX) model and the Eigensystem Realization Algorithm (ERA), and have developed an ARX/ ERA method to identify the coefficient matrices of state equation of motion for structures from measured I/O data.
An ARX model is established from the measured I/O data to obtain the Markov parameters, which is used for identification of coefficient matrices of state equation of motion by ERA.
In the second part, the inverse eigensensitivity method is employed to update the finite element model based on the identified modal parameters. A six degrees of freedom lumped mass system is employed to illustrate the proposed procedure for application in finite element model update.
Keywords: system identification, Auto-regressive with Exogeneous (ARX) Model, Eigen-system Realization Algorithm, inverse eigensensitivity method, model
updating 二、緣由與目的 在結構系統鑑定之領域中， 依訊號處理之域別可 區分為時域與頻域，依待處理之試驗資料又可區分 為(甲)時域之自由響應資料(乙)時域之一般性輸出 入資料(丙)頻域之頻率響應函數資料(丁)時域之脈 衝響應函數(IRF)資料或又稱為馬可夫參數。計畫主 持人所領導之研究小組在前數年之研究工作中已 陸續完成處理甲、乙及丙類資料之探討，本計畫則 針對丁類之 IRF 資料之結構系統鑑定作較深入研 究。FEM 在工程上雖已普遍應用，但仍含帶相當程 度的未確定因素，有需藉由試驗結果來修正，以使 結構分析更精緻。本計畫以結構動態量測識別的模 態參數為參考基礎發展特徵靈敏度反式法之 FEM 結構修正方法。 三、結構之狀態運動方程式 一個 p 自由度系統之運動方程式可以下列之二階微 分矩陣方程式表為： (t) (t) (t) (t) Ew Kw u w
M&& + & + = (1) 其中，M、E與K分別為p×p之質量、阻尼與勁 度矩陣， (t)u 與 (t)w 為p×1之外力與位移向量。令 狀態向量為： = (t) (t) (t) w w z & (2) 則式(1)之二階矩陣微分運動方程式可改寫為下列 之一階矩陣狀態運動方程式： (t) (t) (t) A z Bu z& = c + c (3) 並定義輸出方程式為： (t) (t) (t) C z Du y = + (4) 其中： (t)z 為2p×1之狀態向量， (t)u 為ni×1之輸 入向量、n 為力輸入數， (t)i y 為no×1之輸出向量、 o n 為響應輸出數。A 為c 2p×2p之狀態矩陣，B 為c i n p 2 × 之輸入影響矩陣， C 為no×2p之輸出影響
矩陣，D為no×ni之直接傳輸矩陣。令A 、c B 之c 下標之 c 表示該矩陣為連續時間，而 C 、D在連續 及離散時間系統中完全相同，故不另外表示。A 、c c B 、 C 、D定義如下[1]： − − = − − E M K M I 0 Ac 1 1 (5.a) = −1 c M 0 B (5.b) ] [C C M K C C M E C 1 a v 1 a d− − − − = (5.c) 1 aM C D= − (5.d) d C 、C 、v C 為a p×p之對角線矩陣，對角線之值 為 0 或 1，其不為 0 的位置分別表示位移、速度、 加速度量測所在之自由度。若取 t∆ 為量測採樣週 期，則連續時間狀態運動方程式與輸出方程式可被 轉換成下列的離散時間狀態方程式： + = + = + t t t t d t d 1 t u D z C y u B z A z (6) (7) 其中離散時間狀態方程式的狀態矩陣A 與輸入影d 響矩陣B 使用零階保持方法轉換時分別為：d c c d d c d B I]A [A B A A 1 c t 0 A d B e ) t exp( c − ∆ τ − = τ = ∆ =
∫
(8.a) (8.b) 式(7)沿著 t＝0,1, 2, … ,k 依序展開可得： = + = = + = = = =∑
∑
∑
= − = − − = − − k 1 t k 0 t t k k k t k 1 t k k 1 t t k 1 t k 1 0 1 0 1 0 0 0 u Y u D u B CA y u B A z Du u CB y u B z Du y 0 z d d d d d d d Μ (9) 令： D Y0= ，Y CAd Bd 1 k k= − (10) 式(9)中之最後一行利用Y 矩陣表示輸出向量k y 和k 輸入向量u 之間的關係，即為單位瞬時外力引起之k 動 態 反 應 ， 稱 為 衝 擊 響 應 函 數 （ pulse response function），即u 向量可經由k Y 矩陣的轉換得到k yk 向量。矩陣Y 內包含運動系統的動態特性參數，稱k 為系統馬可夫參數或簡稱為 MP，具備有其值不會 因狀態座標轉換而改變的特性。意即，只要結構系 統確定，Y 應為固定值。k 四、自我迴歸外變數模型 在時間序列模型中，將量測之輸入/輸出資料依時間 先後發生次序以迴歸方法所建立之參數化線性組 合的模型，稱為 ARX 模型。利用此模型可求出自 我迴歸(AR)參數矩陣及外變數(X)參數矩陣，再利用 此二參數矩陣求解 MP 參數。對於一n 組輸入與i no 組輸出之動態系統而言，q 階 ARX 模型可將輸出與 輸入關係表為[2]： t q 0 j j t j q 1 i i t i t y u a y =−∑
+∑
+ = − = − θ φ (11) 式中：下標 t 為離散時間指標，y 為t n 組輸出反o 應，u 為t n 組輸入外力；i a 為雜訊向量，模擬此t 動態系統之量測誤差，假設其為平均值等於零之白 色雜訊；φ 為第 i 個 AR 參數矩陣，i θ 為第 j 個 Xj 參數矩陣。上式為n 組響應輸出o n 組力輸入之 qi 階 ARX 模型，簡稱為 ARX( q , n , o n )模型。將i 上式之時間指標自 t = q 展開至資料點數 N 後，合 併φ 與i θ 成為參數矩陣j P後再整理可得下列之線 性矩陣方程式。 a WP y= + (12) 式 中 ： y 為no×(N−q) 之 輸 出 資 料 矩 陣 ，P 為 ) n q ) n n (( no× o+ i + i 之系統參數矩陣，W 為因次為 ) q N ( ) n q ) n n (( o+ i + i × − 之 系 統 輸 出 入 資 料 矩 陣， a 為no×(N−q)之輸出資料矩陣，這些矩陣之 格式分別表為： ] [yq1 yq2 yN y= + + Λ (13.a) = − − + − − − − q N 1 1 N q N 1 q q N 1 2 N 1 q 1 N q u u u u u u y y y y y y W Λ Μ Ο Μ Λ Λ Λ Μ Ο Μ Λ Λ (13.b) ] [−φ1 Λ −φq θ0 Λ θq = P (13.c) ] [aq1 aq2 aN a= + + Λ (13.d) ARX 系統參數矩陣P可由最小平方法估算： y W P= + (14) 其中，W 代表 W 的擬反式(pseudo inverse)。+整理式(9)與式(11)可以得到 MP 參數矩陣與 AR 參 數矩陣及 X 參數矩陣間的關係如下[3]： 0 0=θ Y (15.a)
∑
= − − = i 1 j j i j i i Y Y θ φ ，i=1,2,… ,q (15.b)∑
= − − = q 1 j j i j i Y Y φ ，i=q+1,q+2,… ∞ (15.c) 前 q 項部份的 MP 參數是分別由 AR 及 X 參數矩陣 決定，其餘的 MP 參數矩陣則是由 AR 參數矩陣來 決定。 五、特徵系統識別運算 ERA 利用已知系統馬可夫參數推算等效狀態系統 矩陣。其運算順序為先建立一個零階 Hankel 矩陣 ) 0 ( H ，再使用奇值分解方法將此矩陣分解，並分析 其奇值以決定系統階數。最後利用一階 Hankel 矩陣 ) 1 ( H 建立一個最小階數的等效系統矩陣。當Y 為k i o n n × 之矩陣時，選取α與β分別作為 Hankel 矩陣 中Y 的區塊列數與區塊行數，可建立一由k α×β個 k Y 所組成之αno×βni的 Hankel 矩陣為[1]： = − − + + + − + + + + − + + 2 â á k á k 1 á k â k 2 k 1 k 1 â k 1 k k ) 1 k ( Y Y Y Y Y Y Y Y Y H Λ Μ Ο Μ Μ Λ Λ (16) ) 1 k ( − H 可分解為三部份： â 1 k á ) 1 k ( P A Q H − = d− (17) 其中P 稱為可觀測矩陣，á Q 稱為可控制矩陣。則â 當 k = 1 時： â á ) 0 ( PQ H = (18) 在αno×βni 之 Hankel 矩陣H(0)中，存在著 o o n n ×α α 正交矩陣R與βni×βni正交矩陣 S ，且 在滿足RTR=RRT =I及 I SS S ST = T = 條件，將 ) 0 ( H 分解如下： T ) 0 ( R S H = Σ (19) 式中Ó為αno×βni奇值矩陣。假設 MP 參數所構成 的H(0)矩陣在無雜訊條件下，奇值矩陣可表為： = 0 0 0 Ón Σ (20) 其中，Ón=diag[σ1,Λ,σi,σi+1,Λ,σn]為由非零奇值 所 構 成 的 對 角 矩 陣 ； 且 奇 值 滿 足 遞 減 關 係 ， 0 n 1 i i 2 1≥σ ≥ ≥σ ≥σ ≥ ≥σ ≥ σ Λ + Λ 。在無雜訊干 擾下，n 為系統自由度 p 之兩倍。若有雜訊干擾時， 則其他本為零之奇值將可能不為零之較小數值，且 依訊號雜訊比而有不同程度之大小。故在決定α及 β大小時必須注意使αno及βni大於或等於 2p。將 式(19)改寫成： T n n n ) 0 ( R Ó S H = (21) n R 及S 為n R及 S 前 n 個行所組成之矩陣且滿足 n T n n n T nR I S S R = = 。比較式(18)及式(21)兩式中右 式，運用平衡識別的概念可令： 2 1 n n á R Ó P = ， T n 2 1 n â Ó S Q = (22) 再令： ] [ o o o o n n n T n I 0 0 E = Λ (23.a) ] [ i i i i n n n T n I 0 0 E = Λ (23.b) 綜合上式(23)、(21)、及式(22)，可導出： 2 1 n n T n 2 1 n (1) − − =Ó R H SÓ Ad (24.a) i n T n 2 1 n S E Ó Bd= (24.b) 2 1 n n T noR Ó E C = (24.c) 又由式(10.a)Y
0=
D
可知： 0 Y D= (25) 式(24.a)、式(24.b)、式(24.c)與式(25)分別構成離散 時間等效狀態方程式相對於式(7)之狀態、輸入影 響、輸出影響與直接傳輸矩陣[Ad,Bd,C ,D ]。將 離散時間狀態矩陣A 進行特徵分析，可將d A 表d 為： 1 d Ø Ã Ø A = − (26) Ø 為含有 2p 個特徵向量之特徵向量矩陣， Ã 為含 有 2p 個特徵值之特徵值矩陣。由二階微分方程式 之特性可知，狀態系統矩陣之系統特徵值皆為成雙 成對出現。每一對共軛複根(λi,λ*i)代表了該系統矩 陣的第 i 個動態模式，可將其所支配的第 i 個自然 振頻與阻尼比分別表為：2 * i i * i i 1 2 * i i 2 * i i i 2 * i i * i i 1 2 * i i i 2 cos 4 ) ln( ) ln( 2 cos 4 ) ln( t 1 λ λ λ + λ + λ λ λ λ = ξ λ λ λ + λ + λ λ ∆ = ω − − (27) (28) 六、特徵靈敏度反式法修正 FEM 模型 本文所使用的特徵靈敏度反式法，乃是在 1974 年 率先由 Collins[4]等人所提出，之後再由 Chen[5]等 人修正成現在所看到的矩陣形式，整個修正程序如 下：(一) 將欲修正的分析模型表成子矩陣的組合。 修正的質量和勁度矩陣可表示為分析模型中子矩 陣的組合，如下式：
∑
= = L 1 i i i U aM M (29.a)∑
= = L 1 i i i U bK K (29.b) 其中：L 為修正係數的數目，a 和i b 為待決定的修i 正係數，如第 i 個修正係數沒有誤差，則a 和i b 之i 值應為 1，但如a 和i b >>1（或 <<1）i ，表示此元素 在建構上有問題。M 和i K 為系統矩陣的子矩陣。i 質 量 矩 陣 和 勁 度 矩 陣 的 第 i 個 macro element matrices M 及i K ，可表示為與個別修正係數有關i 之元素質量矩陣的和，如下式：∑
= = ne 1 j e j i M M (30.a)∑
= = ne 1 j e j i K K (30.b)其中：n 為第 i 個 macro element matrice 中所含質e
量元素的數目， e j M 及 e j K 為第 j 個元素的元素質量 和勁度矩陣。(二) 特徵靈敏度。如果我們將 p 表示 為由修正係數所組成之向量如下：
{
a1 a2 aL b1 b2 bL}
p= Λ Λ (31) 修正模型第 r 個特徵值以下列泰勒級數展開：∑∑
∑
= = = + ∆ ∆ ∂ ∂ λ ∂ + ∆ ∂ λ ∂ + λ = λ L 2 1 j L 2 1 k k j k j Ar 2 L 2 1 j j j Ar Ar Xr p p p p p p Λ Λ (32) 忽略高次項，則上式可以下面近似式表示：∑
= ∆ ∂ λ ∂ = λ ∆ = λ − λ 2L 1 j j j Ar r Ar Xr p p (33) 第 r 個特徵向量的誤差量也可以下式表示：∑
= ∆ ∂ φ ∂ = φ ∆ 2L 1 j j j r A r p p (34) 將式(33)及式( 34)兩式組合成矩陣型式如下： ∆ ∆ ∆ ∆ ∂ φ ∂ ∂ φ ∂ ∂ φ ∂ ∂ φ ∂ ∂ λ ∂ ∂ λ ∂ ∂ λ ∂ ∂ λ ∂ = φ ∆ λ ∆ L 1 L 1 L r A 1 r A 1 r A 1 r A L Ar 1 Ar L Ar 1 Ar r r b b a a b b a a b b a a Μ Μ Λ Λ Λ Λ (35) 或{ }
∆r (n+1)×1=[ ]
Sr (n+1)×2L{ }
∆p 2L×1 (36) 靈敏度矩陣[ ]
S 中的各元素，可由 Fox[6]等人在r 1968 所提出的方法獲得如下： r i T r r r i T r i r p M p K p ∂ φ ∂ φ λ − φ ∂ ∂ φ = ∂ λ ∂ (37)∑
= φ = ∂ φ ∂ N 1 j j i rj i r c p (38) 其中： = φ ∂ ∂ φ − ≠ λ − λ φ ∂ ∂ λ − ∂ ∂ φ = j) (r p M 2 1 j) (r p M p K c r i T j j r r i r i T j i rj 由式(29.a)，我們可得： i i M a M = ∂ ∂ ， 0 b M i = ∂ ∂ (39.a) 0 a K i = ∂ ∂ ， i i K b K = ∂ ∂ (39.b) 將式(39)代入式(37)和式(38)，可得： r i T r r i r M a =−λ φ φ ∂ λ ∂ (40.a) r i T r i r K b =φ φ ∂ λ ∂ (40.b)∑
= φ α = ∂ φ ∂ N 1 j j i rj i r a ， (41.a) = φ φ − ≠ λ − λ φ φ λ − = α j) (r M 2 1 j) (r M r i T j j r r i T j r i rj
∑
= φ β = ∂ φ ∂ N 1 j j i rj i r b ， = ≠ λ − λ φ φ = β j) (r 0 j) (r K j r r i T j i rj (41.b) 如果我們只量測到 m 個模態，式(36)變成：{ }
∆ m(n+1)×1=[ ]
Sm(n+1)×2L{ }
∆p 2L×1 (42) (三) 比例化特徵靈敏度。求解式(42)所會遇到的問 題之一為靈敏度矩陣可能會 ill-conditioned，因為特 徵向量的微分通常遠小於特徵值的微分，可由式(40) 與式(41)看出此特性。 i r r i r p p ∂ φ ∂ λ ≈ ∂ λ ∂ (43) 故可將式(42)改寫如下： { } [ ] ∂ φ ∂ ∂ φ ∂ ∂ φ ∂ ∂ φ ∂ λ ∂ λ ∂ λ ∂ λ ∂ λ ∂ λ ∂ λ ∂ λ ∂ ∂ φ ∂ ∂ φ ∂ ∂ φ ∂ ∂ φ ∂ λ ∂ λ ∂ λ ∂ λ ∂ λ ∂ λ ∂ λ ∂ λ ∂ = φ ∆λ λ ∆ φ ∆λ λ ∆ = ∆ L m A 1 m A L m A 1 m A m L Am m 1 Am m L Am m 1 Am L 1 A 1 1 A L 1 A 1 1 A L 1 1 A 1 1 1 A L L 1 A 1 1 1 A m m m 1 1 1 b b a a b b a a b b a a b b a a ' S ' Λ Λ Λ Λ Μ Μ Μ Μ Μ Μ Μ Μ Μ Μ Μ Μ Λ Λ Λ Λ Μ Μ (44) 即{ }
∆' m(n+1)×1=[ ]
S' m(n+1)×2L{ }
∆p 2L×1 (45) 式(45)為 over-determined 的必須條件為： L 2 ) 1 n ( m + > 或 1 n 2L m + > 七、數值計算例：六質量塊系統 採用一 6 自由度彈簧質塊系統，來對此修正方法進 行模擬，此系統之外觀如圖 1 所示。此系統之質量 矩陣和勁度矩陣分別如下： = 1 0 0 0 0 0 0 2 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 2 M 6 10 2 1 0 0 0 0 1 3 2 0 0 0 0 2 5 1 2 0 0 0 1 2 1 0 0 0 2 1 6 1 0 0 0 0 1 2 K × − − − − − − − − − − − − = 只量測第 1、2、3、4 自由度，即由實驗所得前 4 個模態之特徵值和特徵向量為： = λ 15 . 2389934 0 0 0 0 05 . 1394923 0 0 0 0 11 . 883824 0 0 0 0 15 . 342813 X − − − − − − − = φ 1347 . 0 1928 . 0 0642 . 0 4687 . 0 2983 . 0 7309 . 0 1146 . 0 4692 . 0 0183 . 0 2494 . 0 1921 . 0 3089 . 0 0066 . 0 3158 . 0 8267 . 0 2350 . 0 X 此例中質量和勁度矩陣之子矩陣個數不一致，其子 矩陣分別如下 = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 m M1 ,… , = 6 m 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 M6 = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 k K1 ,… , − − = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 9 k 0 9 k 0 0 0 0 0 0 0 0 0 9 k 0 9 k 0 0 0 0 0 0 0 K9 質 量 與 勁 度 矩 陣 共 有 15 個 參 數 （ m1~m6 ＆ k1~k9），在所有子矩陣都有誤差的情況下進行模 擬，相對於各子矩陣的誤差係數分別為： (0.85, 0.85, 0.9, 0.9, 0.85, 0.9, 0.9, 0.9,… 0.9, 0.9 0.85, 0.9, 0.9, 0.85, 0.85) 修正係數所組成之向量如下：
{
a1,a2,...,a6,b1,b2,...,b9}
此例量測到之模態數 m＝4，如用 4 個模態所構成 之特徵靈敏度矩陣來計算修正係數，在收斂條件訂 為{ }
∆p ≤0.001之下，不論有沒有將特徵靈敏度矩 陣比例化，皆無法收斂。故採用下面途徑來改善收 斂結果，將{ }
∆ 和[ ]
S 補滿至 N 個模態（本例 N＝6）， 使式(42)變成{ }
∆ N(N+1)×1=[ ]
SN(N+1)×2L{ }
∆p 2L×1，其中{ }
∆ 內對應到實驗未量測之自由度則以零取代，則 在相同的收斂條件{ }
∆p ≤0.001下，當特徵靈敏度 矩陣沒有比例化，經過疊代 37 次後收斂，當特徵 靈敏度矩陣有比例化，經過疊代 28 次後收斂，可 見比例化有助於改善收斂速度。 八、結論 本研究整合 ARX 與 ERA 系統鑑定法，能合理有效 的識別出結構系統的特性參數。此法優點可以藉由 ARX 模式直接由量測資料取得馬可夫參數，藉由 ERA 免去一般由 ARX 模型識別模態參數需要降階 的困擾。並將識別結果提供給 IEM 更新運算以修正 六自由度質量塊系統。由數值例結果可知，模態參 數識別的越經精確，其修正後之結果越佳。比例化 特徵靈敏度有助於改善修正過程中之收殮速度。 九、參考文獻1. Juang, J. N., "Applied System Identification", 1994, Prentice Hall.
2. 柯文俊，洪振發，"結構動態系統矩陣鑑定"，中 華民國第六屆振動與噪音工程學會研討會論文 集，pp165-173,八十七年五月。
3. LEE, H. C., M. H. HSAIO, J. K. HUANG and C. W. CHEN 1996 Journal of Vibration and Acoustics 118, pp169-175. Identification of Stochastic System and Controller via Projection Filters. 4. Collins, J. D., G. C. Hart, T. K. Hasselman, and B.
Kennedy, “Statistical Identification of Structures”, AIAA Journal, Vol. 12, No. 2, 1974, pp. 185-190. 5. Chen, J. C. and J. A. Garba, “Analytical Model
Improvement Using Modal Test Results”, AIAA Journal, Vol. 18, No.6, 1980, pp. 684-690.
6. Fox, R. L., and Kapoor, M. P., “Rates of Change of Eigenvalues and Eigenvectors,” AIAA Journal, Vol. 6, No. 12, 1968, pp. 2426-2429. m1=2 m2=2 m3=1 m4=1 m5=2 m6=1 k1=1*106 k2=1*106 k3=1*106 k8=2*106 k9=2*10 6 k4=1*10 6 k 5=2*10 6 k 6=1*10 6 k 7=1*10 6 x1 x2 x3 x4 x5 x6 圖 1. 6 自由度質量塊系統
