希伯特轉換於非線性結構系統之損壞診斷(II)

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

希伯特轉換於非線性結構系統之損壞診斷

計畫類別：整合型計畫

計畫編號：NSC 94－2526－Z－006－016 執行期間：94 年 8 月 1 日至 95 年 10 月 31 日

計畫主持人：洪李陵 共同主持人：

計畫參與人員：黃百誼、黃李暉

成果報告類型(依經費核定清單規定繳交)：完整報告

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

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

處理方式：除產學合作研究計畫、提升產業技術及人才培育研究計 畫、列管計畫及下列情形者外，得立即公開查詢

二年後可公開查詢

執行單位：國立成功大學土木工程學系

中 華 民 國 96 年 1 月 31 日

摘要 摘要 摘要 摘要

線性系統中，有勁度項與阻尼項兩個參數需要識別。以線性時變系統而言，每 一個時間點只有一組含量測反應的運動方程式，故無法識別出勁度項及阻尼項兩組 參數。為了識別這些參數，可對原運動方程組進行希伯特轉換，在每一個時間點上，

再得另一條額外的運動方程式。藉由 Bedrosian 定理，新的運動方程式之勁度項及 阻尼項與原運動方程式相同，故可利用這兩條運動方程式，求解兩個待定系統參數。

本文將介紹運用希伯特轉換系統識別方法，求解線性非時變及線性時變系統之參數 歷時。於模擬系統中，除非參數突然變化，否則識別結果幾乎正確。

關鍵詞：系統識別、希伯特轉換、線性時變系統。

ABSTRACT

There are two parameters, stiffness and damping, to be identified in linear vibration system. For linear time-varying systems, the two parameters can not be identified at each time instant by only using the equation of motion together with measured responses.

In order to identify those parameters, one additional equation of motion is supplied by performing the Hilbert transform to the original equation. With the help of Bedrosian's theorem, the stiffness and damping terms are unchanged at the same time instant and then can be solved by using these two equations simultaneously. This method is shown here to find the parameter time histories of linear time-invariant and time-varying systems.

In simulation study the identified results are almost correct except at the time instants when parameters are changed suddenly.

Key words: System identification, Hilbert transform, Linear time-varying system.

目 目 目 目錄 錄 錄 錄

摘要摘要摘要

摘要...I

ABSTRACT...I

目錄目錄目錄

目錄... II 圖目錄圖目錄圖目錄

圖目錄...III 第第第

第 1 章章章章 前言前言... 1前言前言

1.1

研究動機與目的... 1

1.2

文獻回顧... 1

第第第 第 2 章章章章 希伯特轉換希伯特轉換 ... 2希伯特轉換希伯特轉換

2.1

簡介... 2

2.2

基本定義... 2

2.3

希伯特轉換之解析函數 ... 3

2.4 B EDROSIAN

定理 ... 4

第第第 第 3 章章章章 模擬系統之識別模擬系統之識別... 9模擬系統之識別模擬系統之識別

3.1

單自由度系統之識別 ... 9

3.1.1

線性非時變系統 ... 9

3.1.2

線性時變系統... 11

3.2

多自由度系統之識別 ... 12

第第第 第 4 章章章章 標竿結構之識別標竿結構之識別... 42標竿結構之識別標竿結構之識別

4.1

標竿結構... 42

4.2

標竿結構 D ... 42

4.3

分段參數識別 ... 42

4.4

希伯特系統識別... 43

第第第 第 5 章章章章 結論結論... 56結論結論 參考文獻參考文獻參考文獻 參考文獻 ... 57

附錄附錄附錄 附錄 A 主程式處理流程主程式處理流程... 58主程式處理流程主程式處理流程

A.1

線性內插法(L

INEAR INTERPOLATION OF EXCITATION )... 58

A.2

分段希伯特轉換... 58

圖目錄 圖目錄 圖目錄 圖目錄

圖(2.1) 頻譜不重疊之高低頻訊號...7

圖(2.2) 式(2.21)之有效積分範圍 ...7

圖(2.3) 式(2.22)之有效積分範圍 ...8

圖(3.1) 地表加速度歷時 ...17

圖(3.2) 線性非時變模擬系統之絕對加速度歷時...17

圖(3.3) 線性非時變模擬系統之相對速度歷時...18

圖(3.4) 線性非時變模擬系統之相對位移歷時...18

圖(3.5) 線性非時變模擬系統

^{c t ( )} m

識別值

...19

圖

(3.6)

線性非時變模擬系統

^{k t ( )} m

識別值

...19

圖

(3.7)

相對速度與相對位移為數值積分之

^{c t ( )} m

識別值

...20

圖

(3.8)

相對速度與相對位移為數值積分之

^{k t ( )} m

識別值

...20

圖

(3.9) c t ( ) m

識別值之相對誤差

...21

圖

(3.10) k t ( ) m

識別值之相對誤差

...21

圖

(3.11)

線性時變系統狀況一之

^{c t ( )} m ...22

圖

(3.12)

線性時變系統狀況一之

^{k t ( )} m ...22

圖

(3.13)

線性時變系統狀況二之

^{c t ( )}

m ...23

圖

(3.14)

線性時變系統狀況二之

^{k t ( )}

m ...23

圖

(3.15)

線性時變系統狀況三之

^{c t ( )}

m ...24

圖

(3.16)

線性時變系統狀況三之

^{k t ( )}

m ...24

圖

(3.17)

線性時變系統狀況一之

^{c t ( )}

m

識別值

...25

圖

(3.18)

線性時變系統狀況一之

^{k t ( )}

m

識別值

...25

圖

(3.19)

線性時變系統狀況二之

^{c t ( )}

m

識別值

...26

圖

(3.20)

線性時變系統狀況二之

^{k t ( )}

m

識別值

...26

圖

(3.21)

線性時變系統狀況三之

^{c t ( )}

m

識別值

...27

圖

(3.22)

線性時變系統狀況三之

^{k t ( )}

m

識別值

...27

圖

(3.23)

線性時變系統狀況一強動階段

^{c t ( )}

m

識別值之相對誤差

...28

圖

(3.24)

線性時變系統狀況一強動階段

^{k t ( )}

m

識別值之相對誤差

...28

圖

(3.25)

線性時變系統狀況二強動階段

^{c t ( )}

m

識別值之相對誤差

...29

圖

(3.26)

線性時變系統狀況二強動階段

^{k t ( )}

m

識別值之相對誤差

...29

圖

(3.27)

線性時變系統狀況三強動階段

^{c t ( )}

m

識別值之相對誤差

...30

圖

(3.28)

線性時變系統狀況三強動階段

^{k t ( )}

識別值之相對誤差

...30

圖

(3.29)

線性時變系統狀況一相對速度與相對位移為數值積分之

^{c t ( )}

m

識別值

...31

圖

(3.30)

線性時變系統狀況一相對速度與相對位移為數值積分之

^{k t ( )}

m

識別值

...31

圖

(3.31)

線性時變系統狀況二相對速度與相對位移為數值積分之

^{c t ( )}

m

識別值

...32

圖

(3.32)

線性時變系統狀況二相對速度與相對位移為數值積分之

^{k t ( )}

m

識別值

...32

圖

(3.33)

線性時變系統狀況三相對速度與相對位移為數值積分之

^{c t ( )}

m

識別值

...33

圖

(3.34)

線性時變系統狀況三相對速度與相對位移為數值積分之

^{k t ( )}

m

識別值

...33

圖

(3.35)

雙自由度線性非時變系統

¹

1

( ) c t

m

識別值

...34

圖

(3.36)

雙自由度線性非時變系統

²

1

( ) c t

m

識別值

...34

圖

(3.37)

雙自由度線性非時變系統

¹

1

( ) k t

m

識別值

...35

圖

(3.38)

雙自由度線性非時變系統

²

1

( ) k t

m

識別值

...35

圖

(3.39)

雙自由度線性時變系統狀況一

¹

1

( ) c t

m

識別值

...36

圖

(3.40)

雙自由度線性時變系統狀況一

²

1

( ) c t

m

識別值

...36

圖

(3.41)

雙自由度線性時變系統狀況一

¹

1

( ) k t

m

識別值

...37

圖

(3.42)

雙自由度線性時變系統狀況一

²

1

( ) k t

m

識別值

...37

圖

(3.43)

雙自由度線性時變系統狀況二

¹

1

( ) c t

m

識別值

...38

圖

(3.44)

雙自由度線性時變系統狀況二

²

1

( ) c t

m

識別值

...38

圖

(3.45)

雙自由度線性時變系統狀況二

¹ 1 ( ) k t m

識別值

...39

圖

(3.46)

雙自由度線性時變系統狀況二

² 1 ( ) k t m

識別值

...39

圖

(3.47)

雙自由度線性時變系統狀況三

¹ 1 ( ) c t m

識別值

...40

圖

(3.48)

雙自由度線性時變系統狀況三

² 1 ( ) c t m

識別值

...40

圖

(3.49)

雙自由度線性時變系統狀況三

¹ 1 ( ) k t m

識別值

...41

圖

(3.50)

雙自由度線性時變系統狀況三

² 1 ( ) k t m

識別值

...41

圖

(4.1)

標竿結構

D

之外觀

...43

圖

(4.2)

基底和樓版之質心加速度

...44

圖

(4.3)

數值計算和量測質心位移之誤差

...45

圖

(4.4)

樓層勁度項之識別結果

...46

圖

(4.5)

樓層阻尼項之識別結果

...47

圖

(4.6)

樓層質量比之識別結果

...48

圖

(4.7)

樓層勁度項之識別結果（質量比固定為 1）

...49

圖

(4.8)

樓層阻尼項之識別結果（質量比固定為 1）

...50

圖

(4.9)

樓版加速度誤差之比較

...51

圖

(4.10)

樓層勁度項之識別結果（量測加速度）

...52

圖

(4.11)

樓層阻尼項之識別結果（量測加速度）

...53

圖

(4.12)

樓層勁度項之識別結果（量測加速度）

...54

圖

(4.13)

樓層阻尼項之識別結果（量測加速度和位移）

...55

圖

(A.1)

不同

N

值線性內插法所得之系統反應

...61

圖

(A.2) N=1

與

N=10000

所得系統反應之絕對誤差

...62

圖

(A.3)

不同

N

值之相對加速度經數值積分後所得之系統反應

...63

圖

(A.4)

圖

(A.3)

中相對位移之差

...63

圖

(A.5)

分段希伯特轉換示意圖

...64

圖

(A.6)

未疊代前之識別結果

...65

圖

(A.7)

疊代一次之識別結果

...65

圖

(A.8)

疊代四次之識別結果

...66

圖

(A.9)

疊代十次之識別結果

...66

第 第

第 第1章 章 章 章 前言 前言 前言 前言

1.1 研究動機與目的 研究動機與目的 研究動機與目的 研究動機與目的

台灣位於環太平洋地震帶，每當地震發生時，多少伴隨人員傷亡以及建築物受 損。有鑑於此，為減少地震所帶來的災害，深入了解地震並檢討防災措施勢在必行。

民國八十九年九月廿一日清晨，發生了撼動全台的集集大地震，高達芮氏規模

7.3

的 威力重創了台灣，各地建築物受損及人員傷亡不計其數。而在災後，除了清理崩塌建 築物之外，外觀看似完好的建築物，其結構體仍有內部受損的可能性，為能安心地繼 續使用，應建立一套簡易而可靠的方法，以非破壞性檢測的方式勘驗結構體是否有受 損，以保障使用者之生命財產不受威脅，因此系統識別方法之研究確有其必要性。

簡單來說，系統識別乃是利用量測系統輸入與輸出資料，估算系統特性參數，即 勁度、阻尼比及振態。傅立葉轉換是一般較為簡便的系統識別方法，經由分析系統輸 入與輸出量測資料之頻率內涵，識別出系統參數。然而以傅立葉轉換進行系統識別受 許多因素影響，例如所分析的資料是否具有穩態、週期之性質，以及識別對象是否為 線性系統，然則受損結構體所量測到的離散資料都不符合這些要求，因而增加了傅立 葉頻譜分析之誤差。相對於傅立葉頻譜分析，因轉換所採用的核心函數不同，希伯特 轉換對具有瞬時頻率內涵之特性，因此利用希伯特轉換進行系統識別，可得系統參數 在時間域上之瞬時變化，此一特性使希伯特轉換適用於受損結構體之系統識別。本文 將探討希伯特轉換用於系統識別方法之技巧，以期此一方法能更具有實用價值。

1.2 文獻回顧 文獻回顧 文獻回顧 文獻回顧

系統識別可在時間域或頻率域上進行之，傳統上利用傅立葉轉換系統識別乃在頻 率域上分析訊號之頻譜，然而誠如

1.1

中所述，傅立葉轉換有諸多限制，因此一般以 快速傅立葉轉換處理離散資料，並以濾波器或移動窗等分析技巧解決識別時變系統的 問題，然而無法識別時變系統參數瞬時變化之行為，是傅立葉頻譜系統識別方法難以 克服之瓶頸。

由於希伯特轉換對

[1]

具有相同定義域，故可應用至時間域之系統參數識別，黃 鍔等人

[2]

更進一步以經驗模態分解法

(empirical mode decomposition)

搭配希伯特轉 換，建立希伯特－黃轉換

(Hilbert-Huang transform, HHT)

，適用於非穩態、非線性之 訊號分析。畢德成

[5]

利用

HHT

分析集集地震資料，並比較

HHT

與

FFT

系統識別方 法之優劣，顏嘉德

[6]

以

HHT

觀察

921

地震資料三維希伯特頻譜的頻率內涵變化，王 亮元

[3]

以希伯特轉換搭配

Bedrosian

定理，建立快速可靠的系統識別方法。在進階應 用方面，陳宏南

[7]

與吳政憲

[8]

分別將希伯特頻譜應用於橋樑及高樓之非破壞性檢測。

第 第

第 第2章 章 章 章 希伯特轉換 希伯特轉換 希伯特轉換 希伯特轉換

2.1 簡介 簡介 簡介 簡介

希伯特轉換

(Hilbert Transform)

為一種廣泛應用於數位訊號處理

(DSP)

的演算法，

原輸入訊號經希伯特轉換後，可得一複數解析訊號，其實部為原輸入訊號，虛部為希 伯特轉換後之輸出訊號，輸入訊號與輸出訊號為一希伯特轉換對。相較於將序列訊號 在時間域與頻率域之間互換的傅立葉轉換，希伯特轉換對在同一定義域上，構成解析 訊號後，可觀察局部甚至瞬時頻率變化，此一特性對於時變系統之識別相當具有優勢。

2.2 基本定義 基本定義 基本定義 基本定義

定義單維度實數訊號或函數 ( )

u t 之希伯特轉換如以下積分式 [1]

：

1 ( ) 1 ( )

( ) u u

v t P d P d

t t

η η η η

π η π η

∞ ∞

−∞ −∞

=− =

− −

∫ ∫ ^(2.1)

其反轉換為：

1 ( ) 1 ( )

( ) v v

u t P d P d

t t

η η η η

π η π η

∞ ∞

−∞ −∞

= = −

− −

∫ ∫ ^(2.2)

由 於 式 (2.1) 與(2.2)在

η = t

處有 一 奇 異點 ，是 以 將該 積 分 式改 寫 為柯 西 主 值

(principal value of the integral)之形式，以取極限避點積分的方式處理之：

0

1 ( ) 1 ( ) ( )

( ) lim

^A

A A

u u u

v t P d d d

t t t

ε ε ε

η η η η η η

π η π η η

∞ −

−∞ → −

→∞

 

− −

=

∫

− = 

∫

− +

∫

− 

^(2.3)

而式(2.1)與(2.2)所定義之希伯特轉換可簡化表示為摺積關係：

( ) ( ) * 1 v t u t

π t

= (2.4a)

( ) ( ) * 1 u t v t

π t

= − (2.4b)

若由傅立葉轉換對之架構觀之，可將希伯特轉換對表示如下：

( ) ( )

F

u t ⇔ U ω _(2.5a)

( ) ( )

H

u t

⇔

v t _(2.5b)

對函數

¹

π t

施以傅立葉轉換：

1 e

^{i t}

dt i sgn( ) t

ω

ω

π

∞ −

−∞

= −

∫ ^(2.6)

其中定義函數如下：

1 0 sgn( ) 0 0 1 0

ω

ω ω

ω

+ >

 

=  =

 − <



(2.7)

則分別對式

(2.4a)

與

(2.4b)

施以傅立葉轉換可得以下關係式：

( ) ( ) sgn( ) ( )

F

v t

⇔

V

ω = −

i

ω

U

ω

_(2.8a)

( ) ( ) sgn( ) ( )

F

u t

⇔

U

ω =

i

ω

V

ω

_(2.8b)

式

(2.8a)

與

(2.8b)

表示希伯特轉換對經由傅立葉轉換後之相互關係。

2.3 希 希 希 希伯特轉換之解析函數 伯特轉換之解析函數 伯特轉換之解析函數 伯特轉換之解析函數

假設一希伯特轉換對分別為一複數訊號之實部與虛部，則該複數訊號為解析訊 號，最簡單的例子為尤拉公式

(Euler's formula)

：

[ ]

( ) ( )

( )

( ) cos( ) ( )

1 cos( )

( ) let -

cos ( )

1

1 cos( ) sin( )

cos sin

sin

H

u t t v t

v t d y t

t

y t dy y

y y

t P dy t P dy

y y

t

ω

ωη η η

π η ω π

ω ω

ω ω

π ω

∞

−∞

∞

−∞

∞ ∞

−∞ −∞

= ⇔

=− =

−

− +

=

 

=−  − 

 

=

∫

∫

∫ ∫

cos( ) sin( )

H

t t

ω ω

∴ ⇔

_(2.9a)

同理可得：

sin( ) cos( )

H

t t

ω ⇔ − ω

_(2.9b)

因此尤拉公式之實部與虛部為希伯特轉換對，並為解析函數：

( ) ^{z exp} ( ) ^{i t cos} ( ) ^{t i sin} ( ) ^t

ψ = ω = ω + ω (2.10)

同時此一希伯特轉換對還可驗證，對一訊號作希伯特轉換，相當於對此訊號作

90

度的相位移轉：

ω = ω −π 

( )

cos sin

t t

π

2

ω ω 

− =  − 

(2.11b)

假設一解析函數

^ψ ( ) ^t

^，其實部

^{u t} ( )

^與虛部

^{v t} ( )

為希伯特轉換對，則可表示如下：

( ) ( ) ( ) ( ) ^{t u t iv t A t e}

^iφ

^{( )}

^t

ψ = + =

(2.12)

其中

^{A t} ( )

定義為解析訊號之瞬時振幅(instantaneous amplitude)，

^φ ( ) ^t

^{為瞬時相位}

(instantaneous phase)，可分別由以下公式計算之：

( ) ( ) ( )

( ) ( ) ( )

2 2

tan 1

A t u t v t

t v t u t

φ ⁻

 = +



 =  

  



(2.13)

而 隨 著 時 間 序 列 的 變 化 ， 向 量

^{A t} ( )

^於

^{u − v}

^複數 平 面 上 轉 動 的 瞬 時 角 速 度

(instantaneous angular speed) 定 義 為 解 析 訊 號 的 瞬 時 角 頻 率 (instantaneous angular frequency) ω (t )

_：

1

2 2

( ) ( )

( ) ( ) ( ) ( ) ( )

tan ( ) ( ) ( )

d t t dt

d v t u t v t v t u t dt u t u t v t

ω φ

−

=

= = −

+

(2.14)

由於定義瞬時頻率時，並未考慮訊號的頻寬(bandwidth)，然瞬時頻率應為瞬間最 具代表性的頻率，因此較適用於窄頻訊號之分析，若輸入訊號為寬頻訊號，則可利用 經驗模態分解法[2]，將訊號分解為若干窄頻訊號再處理之。

2.4 Bedrosian 定理 定理 定理 定理

若某一訊號為兩個頻率內涵不重疊訊號之乘積，則對該訊號進行希伯特轉換時，

根據 Bedrosian 定理[1]，可以只對高頻訊號進行希伯特轉換，此一特性對於系統識別 具有決定性的影響，在系統特性參數變化平緩的假設之下，系統參數之變化可視為低 頻訊號，故只需對與系統參數交乘的系統反應進行希伯特轉換，系統特性參數可保留 作為未知數。

考慮一包含兩個高低頻的訊號

f (t )

及

g (t )

之訊號

u (t )

：

)

( ) ( )

( t f t g t

u =

         (2.15) 其中

f (t )

為低頻訊號，

g (t )

為高頻訊號。假設

f (t )

與

g (t )

之傅立葉頻譜以ω 為_c 界，如圖(2.1)所示，兩者之頻譜沒有交錯重疊(overlapping)發生，即：

0 ) ( ω =

F 當 ^ω _> ^ω

_c

(2.16a)

0 ) ( ω =

G 當 ^ω _< ^ω

_c

(2.16b)

其中

F ( ^ω )

代表

f (t )

之傅立葉轉換，

G ( ^ω )

代表

g (t )

之傅立葉轉換，ω 為兩者_c 頻率內涵之分界點。根據摺積定理，

u (t )

的傅立葉轉換可表示如下：

∫ − ^∞ ∞

−

⇐⇒ F u G u du t

g t f

F

) ( ) (

) ( )

( ω

     (2.17)

而經希伯特轉換後之傅立葉轉換為：

∫ − ^∞ ∞

−

−

⇐⇒ i F u G u du t

g t f H

F

) ( ) (

) sgn(

)]

( ) (

[ ω ω (2.18)

只對高頻訊號

g (t )

作希伯特轉換，其傅立葉轉換為：

∫ − ^∞ ∞

−

−

⇐⇒ F u i u G u du t

g H t f

F

)]

( ) sgn(

)[

( )]

( [ )

( ω (2.19)

若能證明式(2.18)與式(2.19)之右側相等，則 Bedrosian 定理得證，如下所示：

[ ] [ ^{u ( t ) H f ( t ) g ( t )} ] ^{f ( t ) H} [ ] ^{g ( t )}

H = = (2.20)

將式(2.18)改寫如下：

i sgn( )

ω ^∞

F (

ω

u G u du ) ( )

−

∫

−∞ −

( ) ( ) 2

( ) ( ) 0 2

( ) ( ) 2 0

( ) ( ) 2

c

c c

c c

c c c

c

c

c

c

i F u G u du

i F u G u du

i F u G u du

i F u G u du

ω ω ω ω ω ω ω ω ω ω

ω ω ω ω

ω ω ω

ω ω ω

ω ω ω

ω ω ω

+

−

−

−

− +

−

− − ≥

 

− − ≤ ≤

=  

 − − ≤ ≤

  − ≤ −



∫

∫

∫

∫

(2.21)

而式

(2.19)

可改寫如下：

( ^{0 0} )

( )[ sgn( ) ( )]

( ) ( ) ( ) ( )

F u i u G u du

i F u G u du F u G u du

ω

ω ω

∞

−∞

∞

−∞

− −

= − − −

∫

∫ ∫

( ) ( )

( ) ( )

c c

i F u G u du i F u G u du

ω ω

ω ω

∞

−

−∞

− −

= 

 −



∫

∫ ^(2.22)

顯然式

(2.21)

與式

(2.22)

之被積函數相同，其有效積分範圍分別繪於圖

(2.2)

與圖

(2.3)

。由於式

(2.22)

之積分範圍可再分割如同式

(2.21)

，因此兩積分式完全相同，故式

(2.20)

成立。至此，

Bedrosian

定理已證明對一含有傅立葉頻譜不重疊之高頻及低頻訊

號乘積的合成訊號進行希伯特轉換時，僅需對高頻訊號作希伯特轉換即可，由於時變 系統參數的頻率內涵遠低於系統反應的頻率內涵，故此一定理有助於時變系統的參數 識別。

ω c

ω c

−

( )

F ω

( )

G ω

ω

圖

(2.1)

頻譜不重疊之高低頻訊號

2 ω _c

−

u c

ω − = ω

u c

ω − = − ω

u ω

ω c

ω c

−

ω c

− ω c

2 ω _c

圖(2.2) 式(2.21)之有效積分範圍

○

¹

○

²

○

³

○

⁴

u c

ω − = ω

u c

ω − = − ω

u ω

ω c

ω c

−

ω c

− ω c

圖(2.3) 式(2.22)之有效積分範圍

○

¹

○

²

第 第

第 第3章 章 章 章 模擬系統之識別 模擬系統之識別 模擬系統之識別 模擬系統之識別

3.1 單自由度系統之識別 單自由度系統之識別 單自由度系統之識別 單自由度系統之識別

考慮一單自由度結構物承受地震力作用，其運動方程式如下：

( ) ( ) ( ) ( ) ( )

_g

( )

mu t

+ c t u t
+ k t u t = − mu t

(3.1)

其中 m 、 ( )

c t 、 ( ) k t 分別為該單自由度系統之質量、阻尼與勁度，由於在時變系

統中，系統勁度與阻尼會隨時間改變，故在此以時間

t

之函數表示之。 ( )

u t

、 ( )

u t

、 ( )

u t

分別為相對於基底之加速度、速度、位移歷時，

u t

_g

( )

為基底或地表加速度歷時。

若將式

(3.1)

等號左右兩邊同時除以系統質量

m

，如下式所示，其中 ( )ζ 與

t

ωn

^{( )} t

分

別為時變系統之阻尼比與自然頻率函數：

( ) 2 ( )

_n

( ) ( )

_n

2 ( ) ( )

_g

( )

u t

+ ζ ω

t t u t

+ω

t u t

= −

u t

^(3.2)

單自由度系統依照系統參數 ( )ζ 與

t

ωn

^{( )} t

之行為可分為以下三類：

1.

線性非時變系統：結構行為可用線性之數學方程式模擬，且系統參數 ( )ζ 與

t

n

( ) t

ω 為常數。

2.

線性時變系統：結構行為可用線性之數學方程式模擬，但系統參數 ( )ζ

t

與

n

( ) t

ω 為時間

t

之函數，具時變性。

3.

非線性系統：結構行為具有非線性遲滯行為，回復力為系統反應之非線性函 數，甚至是隱函數。

然 而 以 數 值 方 法 模 擬 複 雜 的 非 線 性 系 統 時 ， 大 多 以 較 簡 單 的 彈 塑 性 行 為

(elastic-plastic behavior)

或雙線性行為

(bilinear behavior)

模擬之，由於是簡化的分析方

法，因此與真實非線性行為有一定的誤差存在。如果非線性系統改以等值線性時變系 統來識別，不僅可得知動態特性的瞬時變化，亦可依動態參數歷時識別結果，再來推 斷系統的非線性行為，因此本文並不就非線性系統討論識別情形，主要著墨於線性時 變系統上。

3.1.1 線性非時變系統

線性非時變系統線性非時變系統線性非時變系統

線性非時變系統之系統參數 ( )ζ 與

t

ωn

^{( )} t

為常數，故將式

(3.2)

改寫如下：

( ) 2

_n

( )

_n

2 ( )

_g

( )

u t

+ ζω

u t

+ω

u t

= −

u t

^(3.3)

或

( ) 2 ( ) 2 ( ) 0

t

n n

u t

+ ζω u t
+ ω u t = (3.4)

其中

u t

^t

( )

為絕對加速度歷時，

u t

^t

( ) = u t

( ) + u t

_g

( )

。

由於系統參數 ( )ζ 與

t

ωn

^{( )} t

為常數，對式

(3.4)

進行希伯特轉換時，僅需對系統反 應 ( )

u t

與 ( )

u t 做運算，如下式所示，其中 v t

^t

( )

、 ( )

v t

與 ( )

v t 分別為 u t

^t

( )

、 ( )

u t

與 ( )

u t 之

希伯特轉換對：

( ) 2 ( ) 2 ( ) 0

t

n n

v t + ζω v t + ω v t =

(3.5)

式(3.4)與式(3.5)為兩條線性獨立方程式，故可解出兩個未知的系統參數 2ζω 與_n

2

ω

n ：

2

( ) ( ) ( ) ( )

2 ( ) ( ) ( ) ( )

( ) ( ) ( ) ( )

( ) ( ) ( ) ( )

t t

n

t t

n

u t v t u t v t u t v t u t v t

u t v t u t v t u t v t u t v t ζω

ω

 = −

 −

 

 −

 = −

 −

(3.6)

為驗證希伯特轉換系統識別方法於線性非時變系統之適用性，假設一單自由度之 線 性 非 時 變 系 統 ， 其 阻 尼 比ζ =

0.025

、 自 然 頻 率ω =_n

20 rps

， 則 系 統 參 數

( ) 2

_n

1

c t m

= ζω = 、

k t m ^{( )} = ω

n

² = ⁴⁰⁰

。該模擬系統受到某一真實系統在地下室所量 測到的地表加速度

u t

_g

( )

作用，其中

u t

_g

( )

已平移至地下室之形心，如圖

(3.1)

所示。假

設地表加速度

u t

_g

( )

在離散點之間呈線性變化，利用線性內插法

(linear interpolation)

， 可得模擬系統之絕對加速度

u t

^t

( )

、相對速度 ( )

u t

與相對位移 ( )

u t 等系統反應，如圖 (3.2)

至圖

(3.4)

所示。

將模擬系統之絕對加速度

u t

^t

( )

、相對速度 ( )

u t

與相對位移 ( )

u t 進行希伯特轉換，

可得其希伯特轉換對

v t

^t

( )

、 ( )

v t

與 ( )

v t ，代入式 (3.6)

可解出該線性非時變系統之系統 參數識別值，如圖

(3.5)

及圖

(3.6)

所示，圖中顯示識別結果幾乎沒有誤差，其微小的誤 差量可能是來自希伯特轉換的過程中之數值誤差。

以希伯特轉換法識別系統參數，前提是必須有地表加速度

u t

_g

( )

、系統反應 ( )

u t

、

( )

u t

與 ( )

u t 之完整量測值，但一般之監測資料多半僅有 u t

_g

( )

與

u t

^t

( )

兩組加速度歷時紀 錄，此時必須對相對加速度

u t

( ) = u t

^t

( ) − u t

_g

( )

進行數值積分，以得到相對速度 ( )

u t

與 與相對位移 ( )

u t 歷時資料。當 ( ) u t

和 ( )

u t 是由模擬系統之 ( ) u t

數值積分而來，則系統

參數識別結果分別如圖

(3.7)

及圖

(3.8)

所示。

當相對速度與相對位移是由數值積分而得時，識別值與真值之相對誤差如圖

(3.9)

及圖

(3.10)

所示。圖中可發現，使用數值積分所得之相對速度與相對位移進行系統參

數識別時，識別結果之誤差明顯增加，尤其在餘震階段地表加速度振幅較小處，由於 加速度與位移之振幅也較小，使式

(3.6)

中分母項 ( ) ( )

u t v t

−

u t v t ( ) ( )

可能趨近於零，表 示兩條方程式近乎平行，因此儘管數值積分與數值模擬所得之系統反應相差並不甚 大，但系統識別的效果卻有相當大的出入。

3.1.2 線性時變系統

線性時變系統線性時變系統線性時變系統

線性時變系統之運動方程式如式

(3.2)

所示，將地表加速度移至左式後可改寫如 下：

( ) 2 ( ) ( ) ( ) 2 ( ) ( ) 0

t

n n

u t

+ ζ ω t t u t
+ ω t u t = (3.7)

由於系統特性參數 ( )ζ 與

t

ωn

^{( )} t

並非常數，在此假設結構物即使受損時，系統特 性參數變化不若系統反應劇烈，故其頻譜分佈於低頻處，與系統反應之頻譜稍有重 疊。為了進行系統識別，在此仍然只對式

(3.7)

中的系統反應 ( )

u t

與 ( )

u t 進行希伯特轉

換，可得第二條運動方程式：

( ) 2 ( ) ( ) ( ) 2 ( ) ( ) 0

t

n n

v t + ζ ω t t v t + ω t v t =

(3.8)

根 據式

(3.6)

， 同 樣 可 解 出 兩 個 未 知 的 系 統 特 性 參 數

^{c t} ( ) ^{m = 2 ζ} ( ) ( ) ^{t ω}

ⁿ

^t

^與

( )

ⁿ

² ( )

k t m

=ω

t

，但因系統反應與系統特性參數之頻譜有所重疊，故需考慮

Bedrosian

定理適用性所造成的誤差，以下藉三種可能發生的線性時變系統做測試。

圖

(3.1)

所示之地表加速度歷時大約在

15

秒處開始進入強動階段，持續約

7.5

秒

後漸漸進入餘動階段。在強震作用下，結構體可能因地震力而產生裂縫，使結構物整 體的勁度降低，另一方面，隔間牆受損所產生的裂縫也可能使結構物變軟，故整體的 阻尼比會升高，因此假設狀況一之線性時變系統，其系統特性參數

^{c t} ( ) ^m

^與

^{k t} ( ) ^m

^在

初動階段為常數，與 3.1.1 節之線性非時變系統相同，經過強動階段之地震力作用後，

( )

c t m

升高為初動階段的 150%、

^{k t} ( ) ^m

降低至初動階段的 50%，在強動階段期間，

系統特性參數均呈線性變化，系統參數歷時分別如圖(3.11)及圖(3.12)所示。

由於進入餘動階段後，受損之結構體可能伴隨重力壓密效應，使整體勁度回升、

阻尼比稍降，故狀況二考慮進入餘動階段後，

^{c t} ( ) ^m

於 7.5 秒內線性降低為初動階段

的 125%、

^{k t} ( ) ^m

線性回升至初動階段的 75%，其歷時如圖(3.13)及圖(3.14)所示。

然而圖(3.1) 中，強動階段之地表加速度並沒有持續增加之趨勢，故狀況三設定 在進入強動階段後，

^{c t} ( ) ^m

於 5 秒內線性上升至初動階段的 150%、

^{k t} ( ) ^m

^線性降低

至初動階段的 50%，持平 5 秒後因重力壓密效應之故，

^{c t} ( ) ^m

於 5 秒內線性降低為 初動階段的 125%、

^{k t} ( ) ^m

線性回升至初動階段的 75%，其歷時如圖(3.15)及圖(3.16) 所示。

以上三種模擬系統之系統參數識別結果如圖(3.17)至圖(3.22)所示，系統識別之過 程詳見附錄。經過分段希伯特轉換、疊代法與高低頻濾波修正等方式處理後，模擬系 統之系統參數識別之效果尚可接受，惟阻尼項誤差稍大，且在疊代過程中，初動與餘 動階段識別值卻產生不少突然跳動之誤差點，因初動與餘動階段的系統反應較小，使 式(3.6)之分母項趨近於零，因此微小的差量即可能造成識別值相當大的誤差，然而這 些誤差點屬於高頻訊號，在系統特性參數為低頻訊號的前提下，可再利用高頻濾波修 正的方式改善之。

圖(3.23)至圖(3.28)擷取模擬系統在強動階段之系統參數識別值相對誤差，圖中可 以看出，在系統特性參數突然改變之際，其誤差格外顯著，即便經過疊代法改善，轉 折處的誤差依然遠大於其他部分，此即為 Bedrosian 定理適用性所造成的誤差。由於 轉折處屬於較高頻之訊號，因此與系統反應之頻譜重疊的程度較其他部分高，所以

Bedrosian 定理較不適用，造成的誤差自然較高。由於疊代過程中是對系統參數 ^{k t} ( ) ^m

之識別值做高頻濾波修正，再利用修正後之

^{k t} ( ) ^m

帶入式(3.7)，可解得

^{c t} ( ) ^m

^，而

為了滿足運動方程式，在此並不再對識別值做修正。

經過模擬系統的驗證，希伯特轉換系統參數識別方法對於時變系統效果尚可接 受，若系統反應由數值積分而來，則系統識別結果如圖(3.29)至圖(3.34)所示，其效果 顯然不若模擬系統好，尤其在餘動階段之識別效果誤差非常大，這部分的誤差可能來 自線性內插法以及數值積分，詳見附錄 A-1。

由於初動階段與餘動階段之識別值一開始就有相當大誤差存在，經疊代後不見得 能改善之，反而可能愈益嚴重，但是疊代法確實可以有效降低強動階段之誤差，因此 如何取捨端看識別之目的為何，本論文主要探討強動階段之系統參數變化，是以初動 與餘動階段之識別效果較差尚可接受。

3.2 多自由度系統之識別 多自由度系統之識別 多自由度系統之識別 多自由度系統之識別

多自由度系統之系統參數識別可視為單自由度系統的延伸，若每個自由度都能得

到系統反應的完全量測值，則每個自由度都可以列出如同式(3.3a)之運動方程式，並 組合成矩陣形式。考慮一個 n 個自由度的剪力屋(shear building)系統，其運動方程式 如下：

u

g

Mu + Cu + Ku = -M1




(3.9)

其中

M

為質量矩陣，

C

為阻尼矩陣、

K

為勁度矩陣，

1

為元素均為 1 之向量，

u

為相對於地表之加速度向量，

u

為相對於地表之速度向量， u 為相對於地表之位移向 量，

_u

_g為地表加速度歷時。

質量矩陣

M

為一對角矩陣，可表示如下：

1

2 0

0

j

n

m m

m

m

 

 

 

 

=  

 

 

 

 

 

M

其中

m

_j代表第

j

個自由度之質量。

若將式

(3.8)

兩邊同乘以質量矩陣之逆矩陣，可將運動方程式改寫如下：

u + Cu + Ku = 0

t

 

^(3.10)

其中

u

為絕對加速度向量，^t

u = u + 1

^t

u

g。

C

矩陣與 K
矩陣為系統參數矩陣，包含元素如下：

1 2 2

1 1

2 3 3

2

2 2 2

1 1

0 0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0 0

j j j j

j j j

n n

n n

c c c

m m

c c c

c

m m m

c c c c

m m m

c c

m m

+ +

+ −

 

 

 

+ −

 − 

 

 

 

 − + − 

 

 

 

 

 − 

 

 

 

C = ^

(3.11)

1 2 2

1 1

2 3 3

2

2 2 2

1 1

0 0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0 0

j j j j

j j j

n n

n n

k k k

m m

k k k

k

m m m

k k k k

m m m

k k

m m

+ +

+ −

 

 

 

+ −

 − 

 

 

 

 − + − 

 

 

 

 

 − 

 

 

 

K = ^

(3.12)

其中

m 、

_j

c 與

_j

k 分別為第

_j

j

個自由度之質量、阻尼與勁度。

若將式

(3.11)

與式

(3.12)

中第一個自由度之質量

m 提出，並假設第 ₁ 1

個自由度與第

j

個自由度之質量比

_{µ = 1}

_j

m m ₁

_j ，系統參數

c

_j

' ₌ c

_j

m ₁

、

k

_j

' ₌ k

_j

m ₁

，則運動方程式

(3.10)

可改寫如下：

ˆ ˆ

u + Cu + Ku = 0

t

^(3.13)

其中 ˆ

C 矩陣與 ˆ K 矩陣之元素如下所示：

1 2 2

12 2 12 2 12 3 12 3

1 1 1 1 1 11

1 1

' ' ' 0 0 0 0

' ' ' ' 0 0 0

0 0 0

ˆ

0 0 ' ' ' ' 0

0 0 0

0 0 0 0 ' '

j j j j j j j j

n n n n

c c c

c c c c

c c c c

c c

µ µ µ µ

µ µ µ µ

µ µ

+ +

+ −

 

 − + − 

 

 

 − + − 

 

 

 

 − 

 

C =

1 2 2

12 2 12 2 12 3 12 3

1 1 1 1 1 1

1 1

' ' ' 0 0 0 0

' ' ' ' 0 0 0

0 0 0

ˆ

0 0 ' ' ' ' 0

0 0 0

0 0 0 0 ' '

j j j j j j j j

n n n n

k k k

k k k k

k k k k

k k

µ µ µ µ

µ µ µ µ

µ µ

+ +

+ −

 

 − + − 

 

 

 − + − 

 

 

 

 − 

 

K =

若能完全量測所有自由度之系統反應，並估計各個自由度之質量比，則由式

(3.13)

可得

n

條方程式，惟尚不足以求解

c 、 ₁ ' c 、…、 ₂ ' c 以及

_n

' k 、 ₁ ' k 、…、 ₂ ' k 等

_n

' 2n

個 未知數，因此與單自由度之系統參數識別方法相同，對式

(3.13)

進行希伯特轉換，可 得另外

n

條方程式：