論文構架

本論文共分為六章，各章內容如下所述：第一章為緒論，介紹研 究動機與目的、文獻回顧及論文架構。第二章將介紹SRIM 識別法，

介紹如何以輸入輸出資料之相關性重建結構系統之狀態系統參數，並 萃取頻率、模態與阻尼比等模態參數的方法與驗證。第三章將介紹狀 態空間 DLV 損傷識別法之原理與應用。第四章以數值模擬分析驗證 狀態空間 DLV 損傷識別法之可行性，考慮結構不同破壞之型式，不 同破壞程度之靈敏度分析。此外，亦針對結構多重破壞以及觀測不足

7

的條件下作分析。第五章以振動台試驗驗證狀態空間 DLV 損傷探測 法於實際應用之可行性。第六章為結論與建議。

8

第二章 SRIM 系統識別分析討論

2.1 前言

無論是結構損傷評估或預估結構受地震作用之行為，都須先求得 結構系統之參數，包括模態參數。吾人可利用系統識別之技巧，經由 實際量測到之結構動態反應與輸入擾動，推算出結構系統的參數，如 自然頻率、阻尼比及模態等，或更進一步萃取出結構之阻尼、勁度矩 陣，作為後續分析之依據。

Juang【20】於 1997 年提出信息矩陣之系統辨識理論(system realization using information matrix；簡稱 SRIM)。此一方法利用資料 之相關性(data correlation)，由輸出與輸入資料在狀態空間系統之架構 下，由可觀測矩陣(observability matrix)與 Toeplitz 矩陣決定結構系統 之狀態空間矩陣(A、B、C 與 D)，可進而推算系統之模態參數，並用 於後續之結構損傷探測分析。

本章將說明如何利用震測資料以 SRIM 法進行系統識別求得結 構系統之A與C矩陣，進而萃取出系統頻率、阻尼比與模態等模態參 數。此外，並探討SRIM 識別時最佳的參數設定以及雜訊之影響。

9

10

r n c ∈R^{2 ×}

B 為連續時間之輸入影響矩陣。

定義輸出向量y(t)，由位移向量x(t)、速度向量x&(t)與加速度向量 )

x&&(t 所組成如下：

y(t)=C_Dx(t)+C_Vx&(t)+C_Ax&&(t) (2.5) 其中，C_D、C_V與C_A分別為R^m×n之位移、速度與加速度輸出影響矩陣。

由式(2.1)求解x&&(t)並代入式(2.5)可得：

[ ]

^{( )}

) (

) ) (

( ^{1 1 1} t

t

t _{D A V A} t C_AM Eu

x Ξ x M C C K M C C

y ^{− −} ⎥ + ⁻

⎦

⎢ ⎤

⎣

− ⎡

−

= & (2.6)

可改寫為：

y(t) =Cz(t)+Du(t) (2.7) 其中，C=[C_D −C_AM⁻¹K C_V −C_AM⁻¹Ξ]，D=C_AM⁻¹E；

) 1

(t ∈R^m×

y 為輸出向量，m 為感應器之數量；

n

Rm^×2

∈

C 為狀態輸出影響矩陣；

r

Rm^×

∈

D 為直接傳輸矩陣。

由於實際量測資料為離散之型式，故須將連續時間之狀態空間方 程式推展為離散時間之型式。令k =kΔt(Δt為取樣週期)，其離散時間 之狀態空間方程式可表示成：

z(k+1) =Az(k)+Bu(k) (2.8) y(k)=Cz(k)+Du(k) (2.9)

11

12

O 為可觀測性矩陣(observability matrix)；

rp

13

N 之最大值為資料長度減 p，Y_p(k)與U_p(k)皆由已知的輸出與輸 入量測資料組成，藉由兩者之相關性可求得O_p矩陣。

各自相關(autocorrelation)與互相關(cross-correlation)矩陣定義如 下：

R_yy =(1/N)Y_p(k)Y_p^T(k) R_yu =(1/N)Y_p(k)U^T_p(k) R_uu =(1/N)U_p(k)U^T_p(k) R_yz =(1/N)Y_p(k)Z^T(k) (2.15) R_zz =(1/N)Z(k)Z^T(k) R_zu =(1/N)Z(k)U^T_p(k) 其中，對稱矩陣R_yy ∈R^mp×mp、R_uu ∈R^rp×rp與R_zz ∈R^2n×2n分別為輸出觀測 矩陣Y_p(k)、輸入矩陣U_p(k)及未知狀態矩陣Z(k)的自相關矩陣；矩陣

rp mp

yu∈R ^×

R 、R_yz∈R^mp×2n與R_zu ∈R^2n×rp分別為輸出觀測矩陣Y_p(k)對於 輸入矩陣U_p(k)、輸出觀測矩陣Y_p(k)對於未知狀態矩陣Z(k)矩陣及未 知狀態矩陣Z(k)矩陣對於輸入矩陣U_p(k)矩陣的互相關矩陣。

2.2.1 萃取系統矩陣

z A、C 矩陣

由式(2.14)左右兩邊乘上(1/N)U^T_p(k)可得：

R_yu =O_pR_zu +T_pR_uu (2.16) 若R_uu⁻¹為非奇異矩陣，則由式(2.16)可得：

T_p =[R_yu −O_pR_zu]R⁻_uu¹ (2.17)

14

同樣地，於式(2.14)左右兩邊乘上(1/N)Y^T_p(k)可得：

R_yy =O_pR^T_yz +T_pR^T_yu (2.18) 又於式(2.14)左右兩邊乘上(1/N)Z^T(k)可得：

R_yz =O_pR_zz +T_pR^T_zu (2.19) 將式(2.17)之T_p代入式(2.18)與式(2.19)移項整理後可得：

R_yy −R_yuR_uu⁻¹R^T_yu =O_p(R_zz −R_zuR_uu⁻¹R^T_zu)O^T_p (2.20) 茲定義R_hh =R_yy −R_yuR_uu⁻¹R^T_yu，R~_zz =R_zz −R_zuR_uu−¹R^T_zu

則式(2.20)可簡化為：

R_hh =O_pR~_zzO^T_p (2.21) 對R_hh∈R^mp×mp作奇異值分解(singular value decomposition；簡稱 SVD)：

[ ] [

n

]

n n n n

n

hh V V U S V

S U S

U USV

R ^T ⎥ ^T =

⎦

⎢ ⎤

⎣

= ⎡

= ₀

0

0 0

0 (2.22)

其中，U_n∈R^mp×2n為R_hhR^T_hh之非零特徵值所對應之左側單位特徵向量；

0

0

n

Rmp^×

∈

U 為R_hhR^T_hh之零特徵值所對應之左側單位特徵向量，

n mp

n₀ = −2 ；

n n n ∈R^{2 ×2}

S 為R^T_hhR_hh之非零奇異值所組成之對角矩陣；

0 0

0

n

Rn ^×

∈

S 為 R^T_hhR_hh 之 理 想 零 奇 異 值 所 組 成 之 對 角 矩 陣 ；

n mp n ∈R ^×2

V 為R^T_hhR_hh之非零特徵值所對應之右側單位特徵向量；

0

0

n

Rmp^×

∈

V 為R^T_hhR_hh之零特徵值所對應之右側單位特徵向量。

15

16

其中，O^∗_p =[O_p(1:(p−1)m,:)^TO_p(1:(p−1)m,:)]⁻¹O_p(1:(p−1)m,:)^T∈R^2n×(p−1)m為

Op 的 偽 逆(pseudoinverse) 矩 陣 ， 而 整 數 p 之 最 小 值 ， 則 需 滿 足 :)

, ) 1 ( : 1

( p m

p −

O 的秩(rank)等於或大於 2n：

2 1

2 ) 1

( − ≥ ⇒ ≥ +

m p n n m

p (2.28) 觀察O_p矩陣可知，其前m 列為即 C 矩陣：

C=O_p(1:m,:) (2.29)

z B、D 矩陣

觀察式(2.13)之T_p矩陣，可知系統矩陣B與D即隱含於其中。

將式(2.16)前乘U^T₀可得：

^{UT0Ryu =UT0OpRzu +UT0TpRuu} (2.30) 由於O_p =U_n，且根據U^T₀ 與U_n之正交性可得：

^{UT0Ryu =UT0TpRuu} (2.31) 將式(2.31)後乘R_uu⁻¹可改寫成：

^{U0TTp =U0TRyuRuu−1} (2.32) 將未知之T_p矩陣分成p 個子矩陣：

^{Tp =[Tp(:,1:r) Tp(:,r+1:2r) L Tp(:,(p−1)r+1:pr)]} (2.33) 其中，由式(2.13)之T_p矩陣與式(2. 23)可推得：

17

18

19

對

A 矩陣進行特徵分析可得：

AΨ=ΨΛ (2.41) 其中，

] [ψ₁,ψ₂,Kψ_2n

=

Ψ ，

⎥⎥

⎥

⎦

⎤

⎢⎢

⎢

⎣

⎡

=

n 2 1

0

0 λ λ

L M M M M M

K

Λ (2.42)

n

R²n^×2

Ψ∈ 為特徵向量組成之矩陣；

n

R²n^×2

Λ∈ 為特徵值所組成之對角矩陣。

特徵向量Ψ為系統之模態向量矩陣，可經由

C 矩陣將其轉換至輸

出座標上，即可求得降階系統之模態向量矩陣：

Φ=CΨ (2.43) 其中，Φ∈R^m×2n。

結構系統之自然頻率與阻尼比，可由Λ_c的實部與虛部求得，將Λ 轉換為Λ_c如下：

c = ln(ΛΔt )

Λ ，Λ_c =diag(λ_c,1,λ_c,2,_Lλ_c,2n) (2.44) 其中，λ_c,i =α_i ± jβ_i =−ξ_iω_i ± jω_i 1−ξ_i² (2.45)

ω 為系統第 i 模態之自然頻率； i

ξ 為系統第 i 模態之阻尼比。 i

由式(2.45)可解得ω_i與ξ_i：

ω_i = α_i² +β_i² (2.46)

20

2 2

i i

i

i α β

ξ α

− +

= (2.47)

由於特徵值Λ與特徵向量Ψ均為共軛複數的形式，因此所求得之 等效自然頻率、阻尼比與模態數量均為系統自由度的兩倍，且以兩兩 共軛成對出現，實際上求得之模態參數仍與自由度的數量相同。

2.3 SRIM 之參數研究

在 SRIM 識別法的所有設定參數中，以 p 與 N 對於計算時間影 響之最直接，因此若能選擇最適當的參數值進行識別分析，則可提升 效率。

由式(2.16)轉換至式(2.17)之必要條件為R_uu須為非奇異矩陣，以 確保其逆矩陣存在存在，因此N 與 p 之選擇須滿足Rank R( _uu)等於 rp。

式(2.28)中，整數 p 之最小值須滿足O_p(1:(p−1)m,:)的秩(rank)大 於或等於2n；亦即：

2 1

2 ) 1

( − ≥ ⇒ ≥ +

m p n n m p

此外，為了瞭解 SRIM 參數的強健性(robustness)，本節將振動反 應之訊號中加入不同噪訊比之雜訊【21】，並考慮不同之 p 值及 N 值進行識別，由模態正交性之良窳判斷識別結果是否精準。本節以一 座五層樓剪力屋架為例進行分析，如圖 2.1 所示。其模態參數由 ETABS 特徵分析得到結果歸納於表 2.1 及 2.2，結構受地震擾動之地

21

震反應歷時如圖2.2 所示。

於原始訊號中分別加入相同比例之噪音，定義二者之比例關係 (noise-to-singal ratio，簡稱NSR)為：

100% RMS

NSR = RMS ×

S

N (2.48)

其中，RMS_S 與RMS_N 分別表示原始訊號及對應之雜訊之均方根值 (root mean square，簡稱RMS)。

因此，可將式(2.8)、(2.9)改寫成含有噪音之訊號：

) ( ) ( ) ( )

1

(k Az k Bu k w k

z + = + + (2.49a) y_v(k)=Cz(k)+Du(k)+v(k) (2.49b) 或

) ( ) ( )

(k u k w k

u_v = + (2.50a) y_v(k) =y(k)+v(k) (2.50b) 其中，y_v(k)∈R^m×1為加入噪音之輸出訊號；

) 1

(k ∈R^m×

v 原各輸出反應對應之噪音；

) 1

(k ∈R^r×

uv 為加入噪音之輸入訊號；

) 1

(k ∈R^r×

w 原各輸入反應對應之噪音；

本文考慮白噪音(white noise)，以MATLAB^®指令randn建立之。

茲考慮地震反應歷時資料受不同程度的雜訊汙染，包括：

22

Case NSR00-無雜訊；Case NSR05-噪訊比 5%；Case NSR10-噪訊比 10%；Case NSR15-噪訊比 15%及 Case NSR20-噪訊比 20%。各案例 之加速度歷時歸納於圖2.2~2.6。

z 固定 N 值，改變 p 值

茲令 N=3000，考慮包括 p=3、4、5、10、30、50、75、100、125 及150 等，針對五種不同雜訊比之分析結果(如圖 2.7~2.11)討論如下：

Case NSR00：在無噪訊汙染之情況下，p=3(最小之 p 值)時即收斂，

可精準識別出五個振態。

Case NSR05：在噪訊比為 5%之情況下，p=3(最小之 p 值)時即收斂，

可精準識別出五個振態。

Case NSR10：在噪訊比為 10%之情況下，p=10 時達到收斂，可精準

識別出前四振態，而第五振態有些許誤差，但當 p=150 時第五振態之誤差變大許多。

Case NSR15：在噪訊比為 15%之情況下，p=5 僅能精確識別出前四

振態。在噪音的影響下高頻變得不易識別，p 值在 10 以後卻僅能識別出前四振態，無論 p 值考慮多大卻無 法將第五振態正確識別出來。

Case NSR20：在噪訊比為 20%之情況下，p=10 僅能精確識別出前四

振態。在噪音的影響下高頻變得不易識別，p 值在 10

23

以後卻僅能識別出前四振態，無論 p 值考慮多大卻無 法將第五振態正確識別出來。

除了由各振態頻率識別結果判斷 SRIM 法之精確性外，本文將進 一步由振態間之正交性良窳進行檢核。在理想的情況下，結構系統的 正交性應滿足：

=0 Φ

Φ^T_iM _j ，i≠ j (2.51) 其中，Φ_i為第 i 振態之模態向量；

Φj為第j 振態之模態向量；

M 為質量矩陣。

茲定義

j T i

ij =Φ MΦ

ε

(2.52) 以及整體振態正交性誤差為

∑

≠

=

j

i ij

n

2

ˆ 1 ε

ε (2.53)

其中，nˆ =C₂ⁿ，n為系統之模態數。本例之結構系統共五個模態，故

=5

n ，nˆ = C₂⁵ =10。

分析結果歸納於表 2.3 與圖 2.12。由於噪訊比達 15%及 20%時，

能正確識別出的正確頻率只有前四組，因此誤差比起噪訊比5%、10%

較大許多。另外，若配合NSR=15%及 20%兩案例，只取前四組有效

24

模態之正交性誤差(如表 2.4 與圖 2.13 所示)檢核其收斂性，可發現當 5

~ 3

p= 時模態正交性較差；當p=10以上時其識別結果已趨穩定。而 模態之正交性誤差隨著噪訊比的比例增大而變大的趨勢。若 p 超過 50 時，其識別結果未必較佳，因此後續分析時均使用p=10。

z 固定 p 值，改變 N 值

由於 N 值越大所需計算的時間越長，因此若能找到達到收斂結 果之最小 N 值，將有助於後續分析時提高計算效率。本節考慮兩種 不同震波為輸入擾動進行分析，包括：

CASE1：Kobe 地震，PGA=0.1g，Δt =0.005sec CASE2：White Noise 地震，PGA=0.1g，Δt =0.01sec CASE3：White Noise 地震，PGA=0.1g，Δt =0.005sec

三種案例均以p=10進行分析，考慮之 N 值包括 N=100、200、

300、400、500、600、700、800、900、1000、1100、1200、1300、

1400、1500、1600、1700、1800、1900、2000、2500、3000、3500、

4000、4500、5000、5500、6000、6500 至 7000。

分析結果歸納於表 2.5~2.7 與圖 2.14。

CASE1：表 2.5 之結果顯示，當 N=100~1500 時，其分析結果誤差較

大，但隨著N 值增加，誤差逐漸縮小。當 N=3000 時已達到 收斂，再增加N 值其識別之結果未必較佳。

25

CASE2：表 2.6 之結果顯示，當 N=100~1000 時，其分析結果誤差較

大，但隨著 N 值增加，誤差逐漸縮小。當 N=1500 時已達到 收斂，再增加N 值其識別之結果未必較佳。

CASE3：表 2.7 之結果顯示，當 N=100~1000 時，其分析結果誤差較

大，但隨著 N 值增加，誤差逐漸縮小。當 N=1500 時已達到 收斂，再增加N 值其識別之結果未必較佳。

由圖 2.14 顯示，整體而言，隨著 N 值增加，誤差逐漸縮小。而 White Noise 之誤差比 Kobe 好，其原因為擾動之能量較平均；由兩種 White Noise 案例中可看出當Δt越小誤差也越小。

由分析結果可發現，基本上只要選擇 N 值=3000 以上即可確定達 到收斂，因此後續分析時將使用N =3000。

26

第三章 DLV 損傷識別分析之理論

3.1 前言

Bernal【12】於 2002 年提出以柔度矩陣之變化為基礎的結構損 傷識別方法，其主要概念為利用破壞前、後結構柔度矩陣之變化，將 其差異矩陣作奇異值分解(singular value decomposition，簡稱 SVD)，

以所得到對應於零特徵值之特徵向量作為荷載，施加於破壞前的結構 上，再由其應力分析結果萃取出最有可能的破壞構件，作為結構損傷 探測之依據。

Bernal【19】並於 2006 年提出新的 DLV 法，在狀態空間系統建 立柔度矩陣，仍以柔度矩陣之變化為基礎，發展出改良式的破壞定位 向量法(本文簡稱狀態空間 DLV 法)。由於狀態系統之參數可直接由 SRIM 識別法出來，並由其中之部分參數建構柔度矩陣，毋須假設質 量矩陣，因此相較於以 ARX 識別結果進行之分析更為精確。Bernal 並依經驗訂定指標作為篩選破壞定位向量及損傷構件之依據，亦針對 平面桁架結構進行分析。結果顯示，狀態空間 DLV 法在多重桿件受 損之情況下可正確定位出來。

27

28

r n c ∈R^{2 ×}

B 為連續時間之輸入影響矩陣。

假設初始條件為零的情況下，且B_c中之E ∈R^n×n 為單位矩陣及 ) 1

(t ∈R^n×

u 的情況下，對式(3.4)作拉普拉斯轉換(Laplace transform)，則 sz(s)=A_cz(s)+B_cu(s) (3.5) 或

z(s)=

(

sI−A_c

)

⁻¹B_cu(s) (3.6)

考慮位移狀態的部分：

x(s) =C₀

(

sI−A_c

)

⁻¹B_cu(s) (3.7) 其中，C₀ =

[

I 0

]

∈R^n×2n為位移狀態輸出矩陣。

當s=0時(靜態)即為結構之柔度矩陣F為，亦即：

F =−C₀A_c⁻¹B_c (3.8)

此外，吾人可由位移輸出方程式推導輸出資料為速度或加速度與

Ac矩陣之關係如下：

x(t) =C₀z(t) (3.9) 將式(3.9)微分並代入式(3.4)可得：

x&(t) =C₀z&(t)=C₀A_cz(t)+C₀B_cu(t) (3.10) 將式(3.10)微分可得：

x&&(t) =C₀A_cz&(t)+C₀B_cu&(t) =C₀A²_cz(t)+C₀A_cB_cu(t)+C₀B_cu&(t) (3.11)

29

其中，

[ ]

0 0

0 ₁

0 ⎥ =

⎦

⎢ ⎤

⎣

= ⎡ ₋

E B M

C _c I (3.12)

因此，式(3.10)與式(3.11)可簡化為：

x&(t) =C₀A_cz(t) (3.13) x&&(t) =C₀A²_cz(t)+C₀A_cB_cu(t) (3.14)

另一方面，x&(t)及x&&(t)亦可分別表示為：

x&(t) =C₁z(t) (3.15) 及

~ ( ) )

( )

(t C₂z t Du t

x&& = + (3.16) 其中，C₁ =

[

0 I

]

∈R^n×2n為速度狀態輸出矩陣；

x&& = + (3.16) 其中，C₁ =

[

0 I

]

∈R^n×2n為速度狀態輸出矩陣；

在文檔中 狀態空間DLV法在剪力構架之地震損傷探測分析與實驗驗證 (頁 29-0)