利用自我迴歸模型與模態曲率理論以研究懸臂樑結構之損傷偵測

(1)

國立台灣大學工學院工程科學及海洋工程學系碩士論文

Graduate Institute of Engineering Science and Ocean Engineering College of Engineering

National Taiwan University Master Thesis

利用自我迴歸模型與模態曲率理論以研究懸臂樑結構之損傷偵測

Using AR Model and Theory of Curvature Modal Shape to Study Damage Detection of Cantilever Beam

鄧皓元 Hao-Yuan Teng

指導教授：柯文俊博士 Advisor: Wen-Jiunn Ko, Ph.D.

中華民國 102 年 11 月

November, 2013

(2)

I

中文摘要

損傷偵測於結構工程問題中為一門重要研究，而樑式結構於工程中又有著廣泛之應用，如：橋梁、建築等。對於樑式結構使用上之安全檢測已存在許多方法；

有鑑於近年來結構系統識別發展的完善，並參考了以結構模態參數作為偵測結構損傷狀況之相關文獻後，其內容多以有限元素法作數值上之模擬，但未有許多考慮實際應用狀況之研究。因此本文將結構系統識別所得之模態參數結合結構系統之模態曲率理論考慮實際狀況下偵測懸臂樑結構損傷之情形。

本文使用自我迴歸模型及狀態空間系統理論實際識別出懸臂樑結構之模態參數，利用上述識別成果估算結構之雷利阻尼矩陣係數，透過此方式建立近似於實際狀況懸臂樑之有限元模型並由電腦程式分別模擬多種受損狀況下受損與未受損懸臂樑結構含雜訊之輸出響應，由此一流程即可模擬在僅含輸出響應之情況下，首先以自我迴歸模型反向識別出懸臂樑之模態參數，並基於結構模態曲率之理論對懸臂樑作損傷狀況之偵測。而在考慮實際情況所假設之各種不同受損狀況中，本文以上述方式模擬之成果顯示此損傷偵測方法於實際情況下對受損懸臂樑結構能給予一定程度之鑑別功效。

關鍵字：系統識別、自我迴歸模型、狀態空間系統、損傷偵測、模態曲率

(3)

II

Abstract

Problems of damage detection in structural engineering are a major research topic. Beams have extensive applications in structures. For example, bridges and buildings are about that. The use of beam structure already exists in many safety testing methods. The development of structural system identification plays a good role in these fields in recent years. With the references of damage detection by structures modal parameters, its content are more about using finite element method to simulate damage detection. This thesis consider a way to detect the injury situation of cantilever structure through modal parameters obtained from system identification and the curvature modal shape theory.

This thesis uses the Autoregressive model and the state-space system to identify modal parameters of cantilever beam structure and estimate its coefficient of Rayleigh damping. With the above steps, this thesis create finite element model to be similar to the actual cantilever beam. By generating responses of damaged and undamaged cantilever beams from a computer program, we can simulate the damage detection based on curvature modal shapes. The results showed that the proposed method for the variety injury situations has a good effect.

Keyword: System Identification, AutoregRessive Model(AR), State-Space System, Damage Detection, Curvature Modal Shape

(4)

III

謝誌

我能順利完成這份研究論文真的需要感謝非常多人在各方面的協助，無論是在專業領域或是陪伴我度過生活上的各種大小雜事等，都是讓我能度過這兩年多研究生涯的動力。

首先承蒙並且感謝各位口試委員：宜蘭大學生物機電工程學系程安邦教授、

澳汰爾工程股份有限公司總經理徐培譽博士、本所陳國在教授、薛文證教授。

感謝您們不吝提出建議與指正讓這份碩士論文更加地完善。

感謝指導教授柯文俊老師，這兩年多來您的諄諄教誨都讓我覺得自己成長不少，您常教導我們嚴謹做人處事的態度我會銘記在心，非常幸運能接受您的指導，除了感謝還是感謝，希望您身體健康順心。

感謝研究室的后鐘學長、偉庭學長及仁豪學長，你們總是在我不知所措時給予協助，感謝你們耐心地傳承實驗室大小事務。感謝旻毅學長，你平時給的種種意見都讓我有不同的思考方向。感謝宜庭及俊蔚，你們的協助使我的口試能順利完成，祝你們未來研究順利。鴻洲，感謝你的各種幫助，祝你一切順利。

感謝好友莊又全，你的數學能力絕對是我心目中的第一。感謝近四年陪伴我的筱迪，你讓我有信心能順利取得碩士學位。

最後特別感謝我的父母親，謝謝您們給了我一個這樣的環境讓我可以不用擔心家中事務專心讀書而能順利完成學業，萬分感謝，希望你們身體健康愉快。

鄧皓元謹識

2013 年 11 月於新北市永和

(5)

IV

圖目錄

圖 2.1 單自由度系統之 ARX 模型特徵根分布圖 ... 18

圖 2.2 狀態空間系統識別流程圖 ... 24

圖 3.1 直立式懸臂樑結構示意圖 ... 26

圖 3.2 數值方法求解特徵方程式示意圖 ... 29

圖 3.3 直立式懸臂鋼樑實際結構圖 ... 31

圖 3.4 直立式懸臂鋼樑實際結構底座 ... 32

圖 3.5 直立式懸臂樑有限元素模型頂部與底部 ... 32

圖 3.6 懸臂樑有限元素模型邊界條件示意圖 ... 33

圖 3.7 加速規黏貼位置示意圖 ... 35

圖 3.8 量測之加速規 ... 35

圖 3.9 加速規量測各通道之加速度歷時圖 ... 36

圖 3.10 加速規量測各通道之加速度時譜圖 ... 37

圖 3.11 AR 模型識別之頻譜模態疊合圖 ... 38

圖 3.12 AR 模型識別之階數與阻尼比關係圖 ... 38

圖 3.13 由 AR 模型識別、理論值及有限元素分析之模態振形比較圖 ... 39

圖 3.14 AR 模型各階數與雷利阻尼係數關係圖 ... 41

圖 3.15 AR 模型各階數之雷利阻尼係數平均值 ... 42

圖 3.16 AR 模型各階數雷利阻尼係數之標準差 ... 42

(9)

VIII

圖 4.1 懸臂樑撓曲示意圖 ... 45

圖 4.2 未受損直立式懸臂樑有限元素模型 ... 52

圖 4.3 直立式懸臂樑分割單元與節點示意圖 ... 52

圖 4.4 直立式懸臂樑截面幾何關係示意圖 ... 53

圖 4.5 直立式懸臂樑分散受損之有限元素模型 ... 55

圖 4.6 前四模態之模態振形(分散多單元受損) ... 56

圖 4.7 前四模態之模態曲率(分散多單元受損) ... 57

圖 4.8 前四模態之模態曲率差(分散多單元受損) ... 57

圖 4.9 直立式懸臂樑集中受損有限元素模型 ... 58

圖 4.10 直立式懸臂樑前四模態之模態振形(集中多單元受損) ... 59

圖 4.11 前四模態之模態曲率(集中多單元受損) ... 60

圖 4.12 前四模態之模態曲率差(集中多單元受損) ... 60

圖 4.13 第一模態之模態曲率差(同單元不同受損程度) ... 62

圖 4.14 第二模態之模態曲率差(同單元不同受損程度) ... 62

圖 4.15 第三模態之模態曲率差(同單元不同受損程度) ... 63

圖 4.16 第四模態之模態曲率差(同單元不同受損程度) ... 63

圖 4.17 基於模態曲率理論之損傷偵測流程圖 ... 65

圖 5.1 原始訊號與含雜訊 SNR 值 30dB 之時間歷時圖 ... 67

圖 5.2 選取通道恰含受損位置 ... 68

圖 5.3 實際選取通道位置示意圖 ... 68

圖 5.4 集中式選取通道位置示意圖 ... 69

(10)

IX

圖 5.5 前四模態之模態振形(單處損傷 D1) ... 71

圖 5.6 前四模態之模態曲率(單處損傷 D₁) ... 71

圖 5.7 第一模態之模態曲率差(單處損傷 D₁) ... 72

圖 5.8 第二模態之模態曲率差(單處損傷 D₁) ... 72

圖 5.9 第三模態之模態曲率差(單處損傷 D₁) ... 73

圖 5.10 第四模態之模態曲率差(單處損傷 D₁) ... 73

圖 5.11 前四模態之模態振形(單處損傷 D₄) ... 74

圖 5.12 前四模態之模態曲率(單處損傷 D4) ... 74

圖 5.13 等分式選取通道第一模態之模態曲率差(單處損傷 D₄) ... 75

圖 5.14 等分式選取通道第二模態之模態曲率差(單處損傷 D₄) ... 75

圖 5.15 等分式選取通道第三模態之模態曲率差(單處損傷 D₄) ... 76

圖 5.16 等分式選取通道第四模態之模態曲率差(單處損傷 D₄) ... 76

圖 5.17 集中式選取通道前四模態之模態曲率差(單處損傷 D₁) ... 78

圖 5.18 集中式選取通道前四模態之模態曲率差(單處損傷 D₄) ... 78

圖 5.19 前四模態之模態振形(多處損傷狀況一) ... 81

圖 5.20 前四模態之模態曲率(多處損傷狀況一) ... 81

圖 5.21 等分式選取通道第一模態之模態曲率差(多處損傷狀況一) ... 82

圖 5.22 等分式選取通道第二模態之模態曲率差(多處損傷狀況一) ... 82

圖 5.23 等分式選取通道第三模態之模態曲率差(多處損傷狀況一) ... 83

圖 5.24 等分式選取通道第四模態之模態曲率差(多處損傷狀況一) ... 83

圖 5.25 集中式選取通道前四模態之模態曲率差(多處損傷狀況一) ... 84

(11)

X

圖 5.26 前四模態之模態振形(多處損傷狀況二) ... 84

圖 5.27 前四模態之模態曲率(多處損傷狀況二) ... 85

圖 5.28 等分式第一模態之模態曲率差(多處損傷狀況二) ... 85

圖 5.29 等分式第二模態之模態曲率差(多處損傷狀況二) ... 86

圖 5.30 等分式第三模態之模態曲率差(多處損傷狀況二) ... 86

圖 5.31 等分式第四模態之模態曲率差(多處損傷狀況二) ... 87

圖 5.32 集中式第一模態之模態曲率差(多處損傷狀況二) ... 87

圖 5.33 集中式第二模態之模態曲率差(多處損傷狀況二) ... 88

圖 5.34 集中式第三模態之模態曲率差(多處損傷狀況二) ... 88

圖 5.35 集中式第四模態之模態曲率差(多處損傷狀況二) ... 89

圖 5.36 前四模態之模態振形(多處損傷狀況三) ... 90

圖 5.37 前四模態之模態曲率(多處損傷狀況三) ... 90

圖 5.38 等分式第一模態之模態曲率差(多處損傷狀況三) ... 91

圖 5.39 等分式第二模態之模態曲率差(多處損傷狀況三) ... 91

圖 5.40 等分式第三模態之模態曲率差(多處損傷狀況三) ... 92

圖 5.41 等分式第四模態之模態曲率差(多處損傷狀況三) ... 92

圖 5.42 集中式第一模態之模態曲率差(多處損傷狀況三) ... 93

圖 5.43 集中式第二模態之模態曲率差(多處損傷狀況三) ... 93

圖 5.44 集中式第三模態之模態曲率差(多處損傷狀況三) ... 94

圖 5.45 集中式第四模態之模態曲率差(多處損傷狀況三) ... 94

圖 5.46 立方雲線(Cubic Spline)示意圖 ... 97

(12)

XI

圖 5.47 前四模態含數值處理後之模態振形(多處損傷狀況一) ... 98

圖 5.48 前四模態含數值處理後之模態曲率(多處損傷狀況一) ... 99

圖 5.49 第一模態經數值處理之模態曲率差(多處損傷狀況一) ... 99

圖 5.50 第二模態經數值處理之模態曲率差(多處損傷狀況一) ... 100

圖 5.51 第三模態經數值處理之模態曲率差(多處損傷狀況一) ... 100

圖 5.52 第四模態經數值處理之模態曲率差(多處損傷狀況一) ... 101

圖 5.53 前四模態含數值處理後之模態振形(多處損傷狀況二) ... 102

圖 5.54 前四模態含數值處理後之模態曲率(多處損傷狀況二) ... 103

圖 5.55 第一模態經數值處理之模態曲率差(多處損傷狀況二) ... 103

圖 5.56 第二模態經數值處理之模態曲率差(多處損傷狀況二) ... 104

圖 5.57 第三模態經數值處理之模態曲率差(多處損傷狀況二) ... 104

圖 5.58 第四模態經數值處理之模態曲率差(多處損傷狀況二) ... 105

圖 5.59 前四模態含數值處理後之模態振形(多處損傷狀況三) ... 106

圖 5.60 前四模態含數值處理後之模態曲率(多處損傷狀況三) ... 107

圖 5.61 第一模態經數值處理之模態曲率差(多處損傷狀況三) ... 107

圖 5.62 第二模態經數值處理之模態曲率差(多處損傷狀況三) ... 108

圖 5.63 第三模態經數值處理之模態曲率差(多處損傷狀況三) ... 108

圖 5.64 第四模態經數值處理之模態曲率差(多處損傷狀況三) ... 109

(13)

XII

表目錄

表 2.1 轉換矩陣對應不同量測資訊之關係 ... 23

表 3.1 直立式懸臂樑結構幾何與材料參數 ... 25

表 3.2 實體元素與樑元素有限元模型前六模態自然頻率(Hz) ... 33

表 3.3 量測之加速規型號與規格 ... 36

表 3.4 AR(40)模型識別與理論解及有限元素分析之模態參數比較 ... 39

表 3.5 AR(40)識別之自然頻率與阻尼比 ... 43

表 4.1 損傷因子之各項參數數值 ... 53

表 4.2 模擬懸臂樑受損狀況 ... 54

表 4.3 懸臂樑分散多單元受損狀況 ... 55

表 4.4 懸臂樑於分散多單元受損前後之自然頻率(Hz) ... 56

表 4.5 懸臂樑集中多單元受損狀況 ... 58

表 4.6 懸臂樑於集中多單元受損前後自然頻率(Hz) ... 59

表 4.7 懸臂樑同單元不同損傷程度受損狀況 ... 61

表 4.8 損傷偵測方法對於損傷識別指標之比較 ... 64

表 5.1 選取通道位置與受損位置之距離 ... 69

表 5.2 模擬單處受損狀況之受損單元與損傷因子 ... 70

表 5.3 單處受損狀況 AR 模型 40 階識別所得自然頻率(單位：Hz) ... 70

表 5.4 模擬多處受損狀況之受損單元與損傷因子 ... 79

(14)

XIII

表 5.5 多處受損狀況由 AR 模型 40 階識別所得自然頻率(單位：Hz) ... 80

表 5.6 多處損傷狀況一之受損單元與損傷因子 ... 98

表 5.7 多處損傷狀況二之受損單元與損傷因子 ... 102

表 5.8 多處損傷狀況三之受損單元與損傷因子 ... 106

表 6.1 本文所模擬各受損狀況之偵測成效 ... 111

表 6.2 本文所考慮各影響損傷偵測因素之結論 ... 111

(15)

XIV

簡稱術語對照表

AR AutoRegressive 自我迴歸

ARX AutoRegressive with eXogeneous 含外變數的自我迴歸

ARMA AutoRegressive Moving Average 自我迴歸移動平均

ARMAX AutoRegressive Moving Average with 含外變數的自我迴

eXogeneous 歸移動平均

NDT Non-Destructive Testing 非破壞性檢測 FEM Finite Element Method 有限元素法 FEA Finite Element Analysis 有限元素分析 CAE Computer Aided Engineering 電腦輔助工程

MS Modal Shape 模態振形

CMS Curvature Modal Shape 模態曲率

(16)

XV

符號說明

a i AR 模型之參數矩陣

A 離散時間狀態系統矩陣

A c 連續時間狀態系統矩陣

A c 等效連續時間之狀態系統矩陣

A s 直立式懸臂鋼樑結構之截面積

bj ARX 模型之外力參數矩陣

b s 直立式懸臂鋼樑結構之寬度

B 離散時間輸入影響矩陣

B c 連續時間輸入影響矩陣

B c 等效離散時間之輸入影響矩陣

B_ 結構系統之應變矩陣

C c 等效離散時間之輸出影響矩陣

C 離散時間輸出影響矩陣

C a 加速度輸出響應矩陣

(17)

XVI

C c 連續時間輸出影響矩陣

C d 位移輸出響應矩陣

C n 懸臂樑結構模態振形理論解公式之係數

C v 速度輸出響應矩陣

C_ 黏性阻尼係數

d 依照選取通道間隔決定之平移距離

D 離散時間直接傳輸矩陣

D c 連續時間直接傳輸矩陣

D c 等效連續時間之直接傳輸矩陣

damage

D 結構之損傷因子

e q 誤差項構成之誤差矩陣

E 結構系統之阻尼矩陣

E s 直立式懸臂鋼樑結構之楊氏模數

, s total

E 結構有限元素模型各元素對應之楊氏模數矩陣

i

 

G x 立方雲線多項式

I 直立式懸臂鋼樑截面之面積慣性矩

(18)

XVII

K 結構系統之勁度矩陣

l i 懸臂樑有限元素模型之元素長度

L 直立式懸臂鋼樑結構之長度

( , )

m x t x 位置與 t 時間樑於m 方向所受彎矩 _x

M 結構系統之質量矩陣

N 結構系統之型函數矩陣

n a AR 及 ARX 模型之 AR 階數

n b AR 及 ARX 模型之 X 階數

O 懸臂樑變形曲線曲率中心

p i 立方雲線多項式常數項係數

 

P x 拉格朗治多項式

Pq ARX 之參數矩陣

q i 立方雲線多項式一次項係數

( )

q t 廣義座標向量

r i 立方雲線多項式二次項係數

s i 立方雲線多項式三次項係數

t s 直立式懸臂鋼樑結構之厚度

(19)

XVIII

T c 轉換矩陣

T cd 僅量測位移響應之轉換矩陣

T cv 僅量測速度響應之轉換矩陣

T ca 僅量測加速度響應之轉換矩陣

s( )

u t 結構系統外力向量

w n 結構系統第 n 模態之自然頻率

( , )

w x t 懸臂樑於 x 位置 t 時間橫向之理論位移

s( )

x t 結構系統位移向量

d( )

y x 懸臂樑變形曲線方程式

y k ARX 模型表示式

m( )

y t 量測所得之輸出向量

y q 輸出響應構成之矩陣

Y q 輸出響應與輸入外力共構成之輸出入資料矩陣

( )

Y x 結構系統受損前後之模態曲率差

( )

z t 狀態變數向量

( )x

 模態振形函數

 懸臂樑彎曲之轉角

(20)

XIX

 雷利阻尼中質量項之係數

a 雷利阻尼中質量項係數之平均值

n 懸臂樑結構各模態模態振形理論解公式之係數

s 雷利阻尼中質量項係數之標準差

 雷利阻尼中勁度項之係數

a 雷利阻尼中勁度項係數之平均值

n 懸臂樑結構各模態模態振形理論解公式之係數

s 雷利阻尼中勁度項係數之標準差

k 模擬量測與計算時所造成之誤差項

( )x

 結構系統之模態曲率

 直立式懸臂鋼樑結構之質量密度

R 懸臂樑變形曲線曲率半徑

 直立式懸臂鋼樑結構之蒲松比

n 結構系統第 n 模態之阻尼比

 系統之特徵根

 系統之特徵向量

var( )x 原始訊號之變異數

var(xxˆ) 原始訊號與量測訊號之變異數

(21)

XX

[UNDCUR] 未受損結構之模態曲率矩陣

[DCUR] 受損結構之模態曲率矩陣

(22)

1

第一章導論

1.1 研究動機

損傷偵測一直是結構、建築工程中的熱門研究領域，損傷偵測之主要目的：

及時發現結構損壞並掌握受損狀況進而避免損傷情況愈趨嚴重輕則造成財產損失，重則足以威脅生命安全。整體來說，對於識別結構損傷狀況的基本問題而言，

主要又分成四項指標：判斷有無損傷、損傷位置、損傷程度及估算壽命。

而由於樑式結構在工程中有著廣泛的應用，因此對於樑式結構的損傷偵測就顯得十分重要；目前對於損傷偵測的方法可以概括為直接法與間接法兩大類型；

直接法方面，主要是利用先進的儀器設備來直接觀測結構的金相變化，即觀察金屬材料組織(或缺陷)及其變化規律的一種材料物理試驗；而間接法方面，則是間接透過量測結構的物理特性，如電阻、聲速、結構的勁度及振動特性等[1]，以達到對結構之損傷能有一定程度的瞭解。

任何一個結構系統的振動問題可由輸入外力、結構系統及輸出響應此三部分構成，且結構在受到輸入外力影響時，因結構的不同，會有獨一無二的輸出響應，

而系統識別(System Identifications)即為一反向問題，一般可分為以下狀況：(1) 同時含有輸入及輸出響應資訊(I/O Data)、(2)僅含有輸出響應資訊(Output Data Only)，而在僅含輸出響應資訊的情況下又可分為：(i)自由振動(Free Vibration Response)、(ii)強迫振動(Forced Vibration Response)。如何由以上資訊反向識別出結構模態參數，如：自然頻率(Natural Frequency)、阻尼比(Damping Ratio)以及模態振型(Modal Shapes)即為系統識別技術所關切。

(23)

2

有鑑於近年來線性結構體系統識別之完善，吾人參考以結構模態參數作為檢測結構損傷準則之相關文獻[1-8, 31]，其相關研究中主要多以有限元素法模擬損傷偵測，而並未找到許多關於考慮實際狀況下偵測結構損傷之研究，故激起吾人將系統識別成果結合模態曲率損傷偵測理論嘗試模擬於實際情況下偵測結構損傷狀況之動機。本文即是要考慮在僅含有輸出響應的情況下，透過量測懸臂樑的輸出響應，以時間序列模型及狀態空間系統識別理論識別出結構物之模態參數，

由此模態參數基於模態曲率之理論模擬懸臂樑結構之損傷偵測。本文研究目的正是期望以間接性的方式偵測出結構損傷，並探討此流程於實際情況中可能之情況。

1.2 文獻回顧

根據結構動力學的內容，對於任意之結構系統而言，其三大特徵參數分別為質量、勁度與阻尼；每個結構系統都有獨一無二的特徵，若掌握此原則，那麼系統在損傷前後的特徵定會有所改變，而這些改變也必將會反應於結構系統的特徵參數以及模態參數中；因此在損傷偵測的相關研究中，依據結構動力特性的變化來檢測結構的受損情況一直受到廣泛的關注。

以結構系統振動特性其相關內容進行量測進而對結構進行損傷偵測在 1960 年代已有學者進行此一類型的相關研究，1969 年 Lifshitz 與 Rotem 透過振動量測以偵測結構損傷[2]。1978 年 Cawley 與 Adam 也以振動量測進行非破壞性之損傷偵測[3]。Stubbs 和 Osegueda 在 1990 年提出了一套以結構受損前後自然頻率變化的靈敏度作為偵測損傷之方法[4]，其中以有限元素法作為分析途徑。

(24)

3

前述偵測損傷之方法皆為比較受損前後之結構模態參數，然而以結構的動態特性進行量測而作為偵測損傷的方式普遍存在著一個問題：無法在同一時間點上分別量測結構在受損前後之數據；就一般而言，這意味著無法使結構在受損前後以完全相同的邊界條件與環境下進行量測，當必須要在破壞前後進行量測時，先天上就產生了一些隨機的誤差進而影響到量測之正確性。1977 年 Doebling 與 Farrar 的研究指出對於在結構受損前後自然頻率改變量所受到隨機誤差來源之影響相對於其他參數是非常小的[5]；但自然頻率之變化對於結構受損前後之靈敏度是不夠理想的，甚至在實驗操作過程中或是受到環境因素影響所造成之誤差即大於此變化，例如：結構因溫度膨脹。因此若要以自然頻率的改變量對結構進行損傷偵測須有非常精確的量測或是結構受損程度很嚴重才會有明顯之鑑別成果。

1980 年後，除了以自然頻率作為結構損傷偵測的標準外，開始有研究學者將結構模態振型應用在結構損傷偵測上；1984 年 West 以模態振型而沒有搭配有限元素之相關方法進行偵測損傷的相關研究[6]，利用 MAC(Modal Assurance Criteria) 指標來計算兩個不同模態振型之相關性，若兩個模態振型分別為由受損前後之結構所得到，那麼此二模態振型之不同處即很有可能是結構受損位置。但若是在對結構系統不夠瞭解的情況下，量測點為接近某模態下的節點，則量測所得資訊也無法清楚反應結構之受損情況。不同於直接以結構模態振型作為損傷偵測之判斷依據，有學者採用模態振型的二次微分，定義上與曲率相同，故其稱為模態曲率 (Curvature Modal Shape)，Pandey 等人於 1991 年以數值模擬的方式證明模態曲率的變化對於結構受損所造成的影響非常敏感[7]，此研究文獻中將模態振型以中央差分法計算得模態曲率，對於文中模擬的懸臂樑以及簡支樑，結果顯示此法用於損傷偵測中效果良好，若模態曲率於受損前後有明顯改變即表示此結構發生受損，同時此明顯改變之位置即為受損處。1996 年 Stubbs 與 Kim 也在沒有對照組的情況下，測試以模態曲率之方法作為偵測結構損傷準則的適用性[8]。Peng 等

(25)

4

人於 2006 年在理論上證明了模態曲率差法對結構損傷偵測之可行性，並提出了計算結構損傷因子矩陣的公式，最後基於兩種線性假設提出了一種簡易評估結構使用壽命的計算方法[9]。相關文獻中也有提到定義損傷指標(Damage Index)的方法[10]。

在結構系統識別的研究領域中，僅含輸出響應可供系統識別分析的情況下，

差分方程式法的使用頗為廣泛；其中包含自我回歸(AutoRegressive, AR)模型、自我迴歸移動平均(AutoRegressive Moving Average, ARMA)模型以及狀態空間模型等。系統識別技術已發展多年，以下針對在僅含有輸出響應資訊的系統識別下作一簡單之文獻回顧：1983 年 Pandit 等人整合了時間序列模型，應用於動態系統分析，其包含了自我迴歸模型、自我迴歸移動平均模型等[11]。1991 年 Pandit 將狀態空間模型歸類為時間序列模型中之特別形式，並將其應用於結構系統之模態與頻譜分析[12]。Overschee 與 Moor 於 1991 年提出了隨機次空間(stochastic subspace)運算法[13, 14]，可由輸出響應以獲得狀態空間模型，其中的物理意義與理論推導本文隨後的章節也有作一詳細的解釋[15]。1996 至 1997 年間，Lardies 利用狀態空間模型及多變數自我迴歸模型識別出在僅有輸出響應情況下之結構模態參數[16, 17]。Hung 與 Ko 等人利用最小平方法原理估計一個多變數自我迴歸模型的自我迴歸參數矩陣，並將識別所得的模態參數應用於結構之修正及更新的問題[18]。而同時具有輸出響應及輸入外力之系統識別中，Ljung 與 Soderstrom 的書籍完整地將含外變數的自我迴歸(Autoregressive with eXogeneous variable, ARX)模型及含外變數的自我迴歸移動平均(Autoregressive Moving Average with eXogeneous variable, ARMAX)模型作理論及參數估計法的介紹[19, 20]，其中又以含外變數的自我迴歸模型之應用最為廣泛。McGraw 等人指出只要適當地提高 AR 模型的階數，在理論上與含外變數的自我迴歸移動平均模型近似等效[21]。

Soderstrom 等人利用最小平方參數估計法來識別一個連續時間中含外變數的自

(26)

5

我迴歸模型[22]。Hung 與 Ko 等人利用最小平方技術識別含外變數的自我迴歸模型，並應用於結構更新及修正問題。除此之外還結合含外變數的自我迴歸模型與 ERA 自混合位移、速度及加速度輸出響應及輸入外力來鑑定結構系統之矩陣[23, 24]。而狀態空間系統萃取模態參數之理論，其中囊括了：離散時間序列模型、

離散狀態空間系統、等效離散狀態空間系統、等效連續狀態空間系統及狀態空間系統其彼此間的轉換過程；此外，為了刪除非系統本身之特徵根而發展出的凝縮理論，及為了萃取系統特徵參數之矩陣轉換法、馬克夫參數不變原理等[25-28]。

本文即是要以在僅含有輸出響應的情況下以自我迴歸模型與狀態空間系統識別出懸臂樑結構之模態參數，以此參數實際利用 Pandey 等人及 Peng 等人提出之模態曲率以及模態曲率差法並考慮於實際結構偵測損傷時可能之情況模擬結構損傷偵測，最後以數值處理方式增進識別成效。

1.3 論文架構

第一章：導論

本章敘述本文之研究目的與研究動機，並對結構損傷偵測研究近況與本文欲使用之結構系統識別等相關研究作一簡單之文獻回顧，希望透過此章節以概略方式介紹後能帶領讀者迅速地瞭解本文之概要。

第二章：時間序列模型及狀態空間系統識別理論

(27)

6

本章將介紹本文所使用之時間序列模型及使用狀態空間系統識別等相關理論。對時間序列模型如何轉換至離散等效狀態空間系統作推導，判斷是否凝縮後再轉換至連續等效狀態空間系統，最後經矩陣轉換法轉換至狀態空間並萃取系統之模態參數。本章末並有一完整之流程圖方便理解完整過程。

第三章：結構系統識別與有限元素模型

由第二章介紹之系統識別理論實際對結構進行識別，並估算系統之雷利阻尼係數以建立有限元素模型；本章也將針對本文所使用之有限元素軟體及有限元素法作一概略之介紹。

第四章：模態曲率於結構損傷偵測之理論基礎

本章首先完整介紹基於結構模態曲率差之損傷偵測理論，包括連續體模態曲率相關理論基礎及多自由度模態曲率有限元素理論。此外也將模擬結構多種損傷狀況，並就此理論實際數值模擬損傷偵測。

第五章：模擬結構輸出響應之損傷偵測

由第三章實際識別出之結構雷利阻尼係數，分別建立完整與未受損之結構，

並在多種損傷狀況下以第四章所介紹之損傷偵測理論對結構受損情況進行偵測，

且嘗試以數值處理方式增進偵測效果。

第六章：結論與未來展望

綜合本文之損傷偵測成果，作出結論並提出未來展望。

(28)

7

第二章時間序列模型及狀態空間系統識別理論

本章將對本文所使用之系統識別相關理論與實踐理論之識別流程作完整說明，其中更詳細的內容可參考文獻資料[24-28]。狀態空間系統是一種用以描述運動方程式之表示式，其系統矩陣中隱含著結構系統之特徵參數；而一組未知的數據資料，排序方式通常為依照時間先後，此種數據資料稱之為時間序列。若能將此時間序列予以合適之模型來描述，如：自我迴歸(AutoRegressive, AR)模型、含外變數的自我迴歸(AutoRegressive with eXogeneous variable, ARX)模型。再透過對其控制、分析並加以預測，利用所得到之差分方程式搭配相關的轉換公式與狀態空間系統即可準確的識別出結構系統未知的特徵與模態參數。

2.1 狀態空間系統

理論上，任何實際結構系統之動態行為皆為無限多自由度之系統，而狀態空間系統即是利用一偏微分方程式做為其控制方程式。然而在系統識別面對的實際 工程問題中，須考量到實驗器材數量之限制，故通常將結構系統視為 n 個自由度 之系統以近似實際結構系統。若以x t 、_s( ) u t 分別表示結構系統_s( ) n1之位移向量及外力向量，則可將結構系統之運動方程式表示為如下之二階微分方程式：

( )_s ( )_s ( ) _s ( )_s

M x t E x t K x t  u t (2.1.1) 其中，M 、E及K分別表示 n 階之質量、阻尼與勁度矩陣； ( )x t 與_s x t 分_s( ) 別為結構之n1速度及加速度響應。式(2.1.1)中，若質量、阻尼與勁度為已知，

(29)

8

即正向地求解一個矩陣系統方程式，然而對於實際結構體而言，要獲取其特徵參數通常不易，結構系統識別即為此類型之反向問題，藉由實驗儀器，如位移計、

速度計與加速規等，量測結構之輸出響應並考慮輸入外力資料以推算結構之模態參數。若將式(2.1.1)所示之二階矩陣微分方程式以一階矩陣微分方程式表示，以

s( )

x t 及x t 表示位移及速度向量作為_s( ) 2n1之狀態變數向量z t( )，則有一階狀態空間運動方程式：

( ) ( )

( )

s s

z t x t

x t

 

  

  (2.1.2)

1 1 1

( ) 0 ( ) 0

( ) ( ) ( )

s n n n n s n n

s

s s

x t I x t

x t M K M E x t M u t

  

  

      

 

       

    (2.1.3)

考慮在實際量測中，l 個自由度且任一自由度只採用一種儀器之情況下(若實 際結構為 n 自由度，通常 ln)，則式(2.1.1)可表示為：

( ) ( ) ( ) ( )

m a s v s d s

y t C x t C x t C x t (2.1.4) 其中，y t 為量測所得之_m( ) l1輸出響應向量，C 、_a C 及_v C 分別為 l n_d  之加速度、速度與位移輸出響應矩陣(Output Influence Matrix)，且這些矩陣中各元素之數值僅為 0 或 1，其意義表示若值為 1 則顯示該自由度使用了其矩陣代表之量測儀器，而若值為 0 則反之。若將式(2.1.4)同樣以一階矩陣微分方程式表示，則有：

1 1 ( ) 1

( ) ( )

( )

s

m d a v a a s

s

y t C C M K C C M E x t C M u t x t

    

 

     (2.1.5)

(30)

9

綜合式(2.1.2)、式(2.1.3)與式(2.1.5)，則有一階狀態空間運動方程式組，將其以矩陣型式表示：

( ) _c ( ) _c _s( )

z t A z t B u t (2.1.6) ( ) ( ) ( )

m c c s

y t C z t D u t (2.1.7) 其中，A 為 2n 階方陣且表示結構狀態方程式之連續時間狀態系統矩陣，_c B_c 為 2n n 且表示連續時間輸入影響矩陣，C 為_c l2n且表示連續時間輸出影響矩陣，D 為 l n_c  且表示連續時間直接傳輸矩陣。在本文往後的章節中，為了方便閱讀及篇幅的簡潔，皆以 A B C D_c, _c, _c, _c表示式(2.1.6)與式(2.1.7)。本節推導之狀態空間方程式組 A B C D_c, _c, _c, _c ，其中系統矩陣A 蘊含了結構之特徵參數，_c 質量、阻尼及勁度，在本章接續的章節中將介紹所使用之矩陣轉換法，將連續等效狀態空間系統轉換至狀態運動方程式，透過凝縮、萃取等過程，目的為識別出系統之線性模態參數。

2.2 時間序列模型

普遍來說，對於觀察所得之數據序列而言，其通常為按照時間排序，此外亦有按照別種方式排列如空間位置，此處將通稱為時間數列。有別於統計學之研究，

時間數列之研究關心的重點在於對觀察所得數值其前後的相關性，雖然某些問題中其觀察所得數值在統計上是獨立存在的，但考慮大多數的情況下，這些數值其前後是有相關性存在的。本節將介紹近年來廣泛應用於系統識別中之時間序列模型。在此附加說明，對於資料統計相關性多為以自我相關函數表示，也因此時間數列分析的方法同樣為建立在測量自我相關函數之基礎上。

(31)

10

2.2.1 ARX 模型

本節中將介紹如何透過輸出、輸入之資料加以預測或作為本文使用系統識別之過程，而本文所主要使用之時間序列模型為 ARX 模型，其表示式為：

1 0

a b

n n

k i k i j k j k

i j

y a y_ b u _ 

 









 (2.2.1)

其中，a 、_i b_j分別表示 l l 之自我迴歸(AR)參數矩陣與 l n 之外力(X)參數矩陣；n 、_a n 則分別表示 AR 與 X 部分之階數；_b _k用以表示模擬時之量測誤差、

計算誤差、未考慮到的高頻模態及可能存在的非線性動態特性。一般 ARX 模型應用於結構模態參數識別時，常與狀態空間系統作結合應用，多假設 AR 與 X 部分之階數相同以簡化問題。此處假設 AR 與 X 部分階數相同，即n_a n_b q，並依序展開式(2.2.1)如下：

1 0 1 1

y b u e

2 1 1 0 2 1 1 2

y a y b u b u e



1 1

1 0

q q

q i q i j q j q

i j

y a y b u e

 

 











1 1 1 1

1 0

q q

q i q i j q j q

i j

y _ a y _{ } b u _{ } e _

 











1 1

2 2 2 2

1 0

q q

q i q i j q j q

i j

y a y b u e

 

     

 













(32)

11

1 0

N N

N i N i j N i N

i j

y a y _ b u _ e

 









 (2.2.2)

上述展開式中，前q項由於初始資訊不足且並未完整地含括整個 ARX 模型應含的參數，故將其捨棄以避免在估計參數時造成負面的影響。而若將 ARX 之參數矩陣以l(lqn q( 1))之P_q表示，則P_q  a₁ a_q b₀ b₁ b_q，並考慮捨棄前q項後之方程式如下：

1 1

1 2

1 2 1 2

1 2

1 1

1 2

q q N

N q

q q n q q q N

q q N

N q

y y y

y y y P e e e

u u u

 

 



   

 

 



 

 

   

     

 

 

(2.2.3)

可將式(2.2.3)以簡潔形式表示：

p q q q

y P Y e (2.2.4)

其中，y 為輸出響應所構成之_q l(Nq)矩陣，Y 為輸出響應與輸入外力所_q

共同構成之(lqn q( 1)) ( Nq)輸出入資料矩陣，e 為誤差項所構成_q l(Nq) 之誤差矩陣。若經過 Moore-Penrose 擬似反式估算後，可以將 ARX 模型之形式替換為如下之系統識別後所得之差分方程系統表示式：

1 0

ˆ ˆ

q q

k i k i j k j

i j

y a y _ b u _

 









(2.2.5)

(33)

12

由本節介紹之內容，在擁有系統之輸出響應及輸入外力後，建立 ARX 模型之 AR 與 X 參數矩陣，便可將其轉換至離散時間之狀態空間模型，以下將接續介紹相關內容。

2.2.2 ARX 模型轉換至離散狀態空間系統

接續前節之內容，將 ARX 模型轉換至離散狀態空間；首先令一狀態向量，

第一個元素為輸出，第二個元素為輸入；並令第k時間選擇前q個輸出響應所構成之向量y_k__1/_{k q}_ 、前q個輸入外力所構成之向量u_k__1/_{k q}_ 如下：

1/

k k q k

k k q

z y u

 

 

  

  (2.2.6)

且

1 1

1/ , 1/

k k

k k q k k q

k q k q

y u

 

   

 

   

   

   

   

   

(2.2.7)

則式(2.2.6)與式(2.2.7)可改寫為如下之狀態空間系統表示式：

1

k k k

z _  Az Bu (2.2.8)

k k k

y Cz Du (2.2.9)

且

(34)

13

1/ 1

1

1/ 1

k k q k

k k q

z y

u

  



  

 

  

  (2.2.10)

11 12

0 22

A A

 

  

  (2.2.11)

1

2

B B B

   

  (2.2.12)

其中， A 為 ARX 模型轉換至離散時間後之q l(  n) q l( n)狀態系統矩陣，

B 則為q l(  n) n之輸入影響矩陣，式(2.2.9)、式(2.2.11)與式(2.2.12)中，其狀態空間參數矩陣表示如下：

1 2 1

11

ˆ ˆ ˆ ˆ

0 0 0

q q

a a a a

A I

I

  

 

 

 

 

(2.2.13)

1 2

12

ˆ ˆ ˆ

0 0 0

b b bq

A

 

 

  

 

 

(2.2.14)

22

0 0 0 0

0 0 0

A I

I

 

 

 

 

 

(2.2.15)

0

1

ˆ 0 0 b B

  

   

  

 

(2.2.16)

(35)

14

2

0 0 I B

  

   

  

 

(2.2.17)

1 ˆ1 ˆ_q

C  a a  (2.2.18)

2 ˆ1 ˆ

C  b bq (2.2.19)

其中，C為轉換至離散時間後之l q l ( n)輸出影響矩陣， D 則為 l n 之直接傳輸矩陣，由本節所介紹之內容，即可將原建立模型後所得之差分方程式轉換至離散時間狀態空間系統 A B C D, , , 。

2.3 凝縮及連續等效狀態空間系統

系統識別中很重要的一環為獲取系統之特徵根，然而在識別過程中，無可避免地會有諸多誤差產生，如：實驗量測誤差、環境造成之雜訊干擾及計算上的誤差等，為求更良好的精確度，在將 ARX 模型透過轉換後得到離散等效狀態空間方程組時，通常會增加 ARX 模型之建模階數，但相對地則會導致因次過大而產生更多非系統之特徵根，因其常不具規則地出現在不同階數之 ARX 模型中，故之為幽靈根(Ghost Roots)，本節即要介紹去除幽靈根並留下系統特徵根之過程，

稱為凝縮。而在離散等效狀態空間方程組經凝縮後，還須將其轉換至連續時間系統上，此一過程可透過連續與離散系統之間的關係式得以完成，以下將依序介紹凝縮及轉換過程。

(36)

15

由上節中 ARX 模型轉換後之離散狀態空間方程組，此未凝縮之狀態空間方程式為：

1

k k k

z _  Az Bu (2.3.1)

k k k

y Cz Du (2.3.2)

將系統之特徵根與特徵向量分別表示為_s及，則 A 矩陣經特徵分解後：

1

1 0

0

ss sp s ss sp

ps pp p ps pp

A



        

              (2.3.3)

式(2.3.3)中，下標s表示系統(System)本身之特徵根， p表示虛假(Spurious) 之特徵根，將 A 矩陣中，系統本身和虛假的所分別對應之特徵向量與特徵根分開 後，再將其代入式(2.3.1)：

1 1

k k k

z _  ^z Bu (2.3.4)

將^¹由左向右乘上式(2.3.4)，並設z_k  ^¹z_k，則有：

1 1

k k k k k

z _    z ^ Bu   z Bu (2.3.5)

以此類推，式(2.3.2)同樣可改寫為：

k k k k k

y   C z Du Cz Du (2.3.6)

將轉換過之狀態向量與輸出訊號同樣分成系統自身及虛假之特徵根與特徵向量，則式(2.3.5)與式(2.3.6)：

(37)

16

1

0 0

s s

k k

p p k

p p

k k

z z B

B u

z z



  

    

   

       

      (2.3.7)

以及

s k

k s p p k

k

y C C z Du

z

  

  

  (2.3.8)

其中

1 s

p

B B B

  

   

   (2.3.9)

C = C_s C_p (2.3.10)

經過篩選，將不屬於系統自身之特徵根與特徵向量忽略，其篩選方式將稍後介紹，則式(2.3.7)及式(2.3.8)有如下之近似式：

1

s s

k s k s k

z _   z B u (2.3.11)

s

k s k k

y C z Du (2.3.12)

而狀態向量轉換前後關係則為：

s s

ss sp

k k

p p

ps pp

k k

z z

 

    

    

    (2.3.13)

由式(2.3.13)可得z_k^s  _ssz_k^s或z_k^s  ^¹z_k^s之關係式，並將此關係式代入式 (2.3.11)與式(2.3.12)：

(38)

17

1 1

1

s s

sszk s sszk B us k

 

      (2.3.14)

1 s s

k s ss k k k k

y  C ^ z Du Cz Du (2.3.15)

接著將式(2.3.14)自左向右乘上_ss，則有：

1 1

s s s

k ss s ss k ss s k k k

z _    ^ z   B u  Az Bu (2.3.16)

式(2.3.15)與式(2.3.16)即為凝縮後之狀態空間系統：

1

s s

k k k

z _  Az Bu (2.3.17)

s

k k k

y Cz Du (2.3.18)

由前節所介紹之形式，可知凝縮後離散時間等效狀態空間系統A B C D, , ,  之轉換式分別為：

1 ss s ss

A   ^ (2.3.19)

ss s

B  B (2.3.20)

1 s ss

C  C ^ (2.3.21)

DD (2.3.22)

而在得到離散時間等效狀態方程組A B C D, , , 後，須將其轉換至連續時間等效狀態方程組A B C D_c, _c, _c, _c ，其轉換關係式：

 

ln /

Ac  A t (2.3.23)

(39)

18

ln( )( ) 1 /

Bc  A AI ^ B t (2.3.24)

Cc C (2.3.25)

Dc D (2.3.26) 前述篩選系統本身之特徵根之過程，即為本節所介紹之凝縮。此法之重點在於如何排除非系統特徵根並保留自身特徵根；事實上，對於任一系統而言，其特徵有一特性，無論使用模型之階數為何，其特徵根為固定值，圖(2.1)將說明此情況：

圖 2.1 單自由度系統之 ARX 模型特徵根分布圖

-1 -0.5 0 0.5 1

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

Real Part

Imaginary Part

ARX order:5 ARX order:10 ARX order:15

(40)

19

假設一單自由度系統，建立其 ARX 模型以識別其特徵根，依序分別為 5 階、

10 階及 15 階，並將各不同階數之特徵根如圖中註記之符號繪於單位圓內，觀察其分布狀況。由圖中所示，不同階數之特徵根將不固定位置地散布於單位圓內，

唯獨圖中右側兩處之特徵根在階數不同時皆重合，若增加更高階數之根值也將重合，而此二位置即為系統自身之特徵根。我們即可利用系統此一特性，配合前述介紹之理論，篩選掉系統之幽靈根藉以選取出系統之特徵根，以達到凝縮之功效。

2.4 矩陣轉換法

透過前幾節所介紹之內容並依其步驟，已可得到一組連續時間等效狀態方程組A B C D_c, _c, _c, _c ，但若欲將狀態空間系統中之實際特徵參數萃取出，還須進行本節將介紹之矩陣轉換法理論，如此一來即可完整識別出系統之線性特徵參數。

本節中，將說明一系列之推導，如何將此一連續時間等效狀態方程組轉換至狀態空間方程組A B C D_c, _c, _c, _c 。

為便於理論推導，假設一 2n2n之非奇轉換矩陣以T 表示，其可使欲轉換_c 之兩組狀態方程組，A B C D_c, _c, _c, _c 與A B C D_c, _c, _c, _c 具如下之關係：

1

c c c c

A T A T^ (2.4.1)

c c c

B T B (2.4.2)

(41)

20

1

c c c

C C T^ (2.4.3)

Dc D (2.4.4) 將此一轉換矩陣T 列分割為上下兩部分，其分別具有相同之列因次： _c

1

2 c c

c

T T T

   

  (2.4.5)

將T 由右向左乘上式(2.4.1)，則有： _c

c c c c

A T T A (2.4.6)

若以矩陣形式表示則為：

1 1

2 2

0_n _n _c _c

c

c c

T T

I A

T T

M K^ M E^

   

    

  

      (2.4.7)

其中，式(2.4.7)可表示為二方程式：

2 1

c c c

T T A (2.4.8)

1 1 2

1 2 2 1

c c c c c c

M KT^ M ET^ T A T A

    (2.4.9)

由式(2.4.8)與式(2.4.9)可以得到T 與_c₁ T 之關係，並由式(2.4.2)： _c₂

1 1

2

0_n _c

c c

c

B T B

T M^

   

   

    (2.4.10)

同樣將式(2.4.10)表示為兩方程式：

利用自我迴歸模型與模態曲率理論以研究懸臂樑結構之損傷偵測

國立台灣大學工學院工程科學及海洋工程學系 碩士論文

Graduate Institute of Engineering Science and Ocean Engineering College of Engineering

National Taiwan University Master Thesis

利用自我迴歸模型與模態曲率理論 以研究懸臂樑結構之損傷偵測

Using AR Model and Theory of Curvature Modal Shape to Study Damage Detection of Cantilever Beam

鄧皓元 Hao-Yuan Teng

指導教授：柯文俊 博士 Advisor: Wen-Jiunn Ko, Ph.D.

中華民國 102 年 11 月

November, 2013

中文摘要

Abstract

謝誌

目錄

圖目錄

表目錄

簡稱術語對照表

符號說明

 

 

第一章 導論

1.1 研究動機

1.2 文獻回顧

1.3 論文架構

第二章 時間序列模型及狀態空間系統識別理論

2.1 狀態空間系統

2.2 時間序列模型

2.2.1 ARX 模型

























2.2.2 ARX 模型轉換至離散狀態空間系統

2.3 凝縮及連續等效狀態空間系統

 

2.4 矩陣轉換法

國立台灣大學工學院工程科學及海洋工程學系碩士論文

利用自我迴歸模型與模態曲率理論以研究懸臂樑結構之損傷偵測

指導教授：柯文俊博士 Advisor: Wen-Jiunn Ko, Ph.D.

第一章導論

第二章時間序列模型及狀態空間系統識別理論