國 立 交 通 大 學
土 木 工 程 學 系
碩 士 論 文
數種濾波技術於線性系統識別之可行性探討
A Study on the Applicability of Several Filtering Techniques
to the Identification of Linear Systems
研 究 生：楊維莘
指導教授：黃炯憲 博士
中 華 民 國 一百零一 年 一 月
i
數種濾波技術於線性系統識別之可行性探討
A Study on the Applicability of Several Filtering Techniques
to the Identification of Linear Systems
研究生：楊維莘
Student：Wei-Hsin Yang
指導教授：黃炯憲 博士
Advisor：Dr. Chiung-Shiann Huang
國 立 交 通 大 學
土 木 工 程 學 系
碩 士 論 文
A Thesis
Submitted to Department of Civil Engineering
College of Engineering
National Chiao-Tung University
in Partial Fulfillment of the Requirements
for the Degree of
Master of Science
in
Civil Engineering
January 2012
Hsinchu, Taiwan, Republic of China
ii
數種濾波技術於線性系統識別之可行性探討
研究生：楊維莘
指導教授：黃炯憲 博士
國立交通大學土木工程學系
摘 要
本研究旨在探討數種常用濾波技術對於非時變系統與時變系統之系統識別模型 參數的影響。本研究選用三種濾波技術來進行比較，分別是卡爾曼濾波演算法(Kalman filtering algorithm)、離散小波濾波法(discrete wavelet de-noising)，以及中位值濾波器 (median filter)。研究中將數值模擬地震反應分為非時變系統與時變系統兩大部分：非 時變系統中，將數值模擬七個自由度的系統之地震反應，利用連續小波轉換建構 ARX (Autoregressive with exogenous input) 模式，進一步估算其系統模態參數。為探討不同 噪訊含量的訊號中，經由上述各種濾波法處理後的差異性，我們將數值模擬地震反應 加入不同比例的白色噪訊，使其噪訊比分別為 10%、20%、30%、40%，再利用濾波 法對含噪訊之地震反應進行濾波處理，而濾波方法對於提升系統識別之效率的可行性， 乃取決於識別出其準確之模態參數時所需的最小階數。時變系統部分，本研究採用數 值模擬五個自由度的系統受地震反應時的訊號，不同於非時變系統，此系統的勁度與 模態的阻尼比與時間相關，利用移動最小平方差法(Moving Least Square method，MLS) 架構 TVARX 模型( Time Varying Autoregressive with exogenous input)，識別其系統之瞬 時模態特性。對於時變系統，其研究方法與探討非時變系統雷同，將上述的濾波技術 應用於處理時變系統中其噪訊比分別為 1%、3%、5%的地震反應，用來評估其濾波技 術對於改善系統識別之效能的可行性。iii
A Study on the Applicability of Several Filtering
Techniques to the Identification of Linear Systems
Student：Wei-Hsin Yang
Advisor：Dr. Chiung-Shiann Huang
Department of Civil Engineering National Chiao-Tung University
Abstract
The main purpose of this study is to investigate the effects that several popular de-noising techniques on identifying the modal parameters of time invariant and time varying linear systems. There are three de-noising techniques under consideration, and they are Kalman filtering algorithm, discrete wavelet de-noising, median filter. When studying invariant time systems, we numerically simulate acceleration responses of a seven degrees of freedom (DOF) system under earthquake. A continuous wavelet transform combined with ARX (Auto-regressive with exogenous input) model is applied to process the acceleration responses and estimate the modal parameters of the system. White noise with different noise levels (noise-to-signal ratio (NSR) =10%, 20%,30%, 40%) is added to the acceleration responses. The de-noising techniques are employed to the noisy responses. The feasibility of these de-noising techniques on handling the noisy data is determined by the least order of ARX yielding accurate modal parameter identification. When investigating time varying systems, we numerically simulate acceleration responses of a five DOF system under earthquake. Their stiffness and modal damping ratios are time dependent. A moving least square method is adopted to construct TVARX model (time varying auto-regressive with exogenous input) and to identify the instantaneous modal parameters of the five DOF system. Similar to the study for a time invariant system those de-noising techniques are utilized to process the noisy acceleration responses with NSR equal to 1%, 3%, and 5%, and their effectiveness on improving the modal identification is evaluated.
iv
誌謝
首先，最要感激的人就是指導教授 黃炯憲老師對我的耐心指導，更感謝老師一 直以來對我的寬容、體諒及指教，讓我在生活上及課業上都得到極大的支持與幫助。 老師對於學術上的熱誠與對學生的關懷令我敬佩萬分，與老師的氣度相比，實在令我 自己感到汗顏，期許自己能努力向老師看齊，學習這種敬業與認真的態度，對於將來 運用於人生的任何一個面向，必將是極大的助益，對老師的感激之情無法言喻，難以 用筆墨來形容。 特別要感謝 翁正強老師對我的教導與啟發，老師諄諄的教誨，予我生活上的支 持、鼓勵與建議，讓我對人生的方向有了更明確的目標，對自我的定位更加清晰。 感謝口試委員 鄭復平老師、 劉俊秀老師、 林昌佑老師的用心傾聽、指導與完 善的建議，讓這份論文更趨於完整。感謝大學時期的導師 王維志老師，指引我未來的 方向；感謝交通大學以及土木系所給予我完善的教育資源與環境，從大學時期到研究 所，這幾年在交大的生活讓我在心靈和教育上得到良好的滋養，得到前所未有的體驗 與快樂。當然，也要感謝清華大學，慷慨地分享資源，帶給我生活上的便利，期望兩 校發展日益精進。 感謝同窗好友鈞誠、江祥、宣妤、連峰、中原、博然、進順、承哲、自勝、柏霖 及其他研究所同學們；學長姐景裕、祖涵、意晴、柏安、威智、靖俞、穎泰、政甯、 明儒等；學弟妹芳琳、裕鈞、旭進、宣治、俊超等，有了你們的陪伴讓我的研究所生 活多采多姿，在學業和生活上都得到大家的照顧，也得到了極其珍貴的友情，更懂得 如何待人處事與良好地溝通。 感謝我的多年好友凱復、正淇、聖文、兆民、宗佑、聖仁、國豪、韋菱、宜萍、 亭宜、祐甄、吉鴻、強瑞，多虧你們，讓我擁有美好的回憶；感謝我的好姊妹們，雅 婷、晏嘉、岱瀚、康亭及其他高中同學，很高興經過這麼多年我們的情誼仍然如此深 厚，妳們是我的支柱。這些年有大家的真摯的友情與情義相挺，分享生活中的心情與 大小事，使我感覺充滿力量，心智與個性也更加成熟。 最後要感謝一直以來最支持我的父親、母親、弟弟允中及喆銘，還有很幸運地在 我畢業前臨門一腳趕上的君瀚，你我都是幸運的人。感謝所有愛護我的家人，感謝整 個大家族，感謝我們擁有彼此，你們都是最疼愛我的家人，更是我最親密的伴侶，有 了你們的照顧，給予我圓滿的愛，讓我一直在備受寵愛的環境裡成長，使我擁有完整 的幸福，才有完整的我。有了你們，我的人生更加美好。感謝所有人對我的祝福。 這些年遇到了很多人、很多事、很多物，促使我在個性、態度、處事及觀念上都v 學習與成長了很多，我的視野更加開闊，心胸更寬大，大大小小的人、事、物，成就 了今天的我，也成就了我的今天，更成就了我的明天。我真是個幸運兒，要感謝的實 在是太多了，就感謝這個世界吧！我真心相信，懷著感恩的心，每一件事堆砌出來的 路，都是在為我美好的未來鋪設更舒適的道路。感恩這個美好的人生，感謝我自己， 盡情自在地享受這一切最美滿的幸福。
vi
目錄
中文摘要 ... ii 英文摘要 ... iii 誌謝……….. ... iv 目錄…….. ... vi 表目錄….. ... viii 圖目錄….. ... ix 第一章 緒論 ... 1 1.1 前言 ... 1 1.2 文獻回顧 ... 2 1.3 研究內容 ... 3 第二章 濾波方法介紹 ... 4 2.1 前言 ... 4 2.2 卡爾曼濾波演算法之介紹 ... 4 2.2.1 卡氏濾波器的基本動態系統模型 ... 5 2.2.2 卡氏濾波器的推導及運算 ... 5 2.3 中位值濾波法之介紹 ... 7 2.3.1 中位值濾波法之公式推導 ... 8 2.4 離散小波轉換之介紹 ... 9 2.4.1 離散小波變換的推導 ... 9 2.4.2 離散小波在濾波上的應用 ... 11 第三章 架構 ARX 模型以及 TVARX 模型 ... 14 3.1 前言 ... 14 3.2 連續小波變換之理論基礎 ... 14 3.3 建構 ARX 模型 ... 17 3.4 以移動最小平方差法架構 TVARX 模型 ... 19 3.5 瞬時模態參數估算 ... 23 第四章 數值模擬之系統識別結果分析 ... 26 4.1 前言 ... 26 4.2 非時變系統之數值模擬的系統識別分析結果 ... 26 4.2.1 小波函數之影響 ... 26vii 4.2.2 無噪訊反應之識別結果 ... 27 4.2.3 含噪訊反應對於識別結果之影響 ... 27 4.3 濾波後對於非時變系統的數值模擬反應識別結果之影響與比較 ... 28 4.3.1 卡氏濾波器濾波後的識別結果 ... 28 4.3.2 中值濾波器濾波後的識別結果 ... 29 4.3.3 離散小波濾波後的識別結果 ... 30 4.3.4 濾波後對於非時變系統的數值模擬反應識別分析的綜合比較 ... 31 4.4 時變系統之數值模擬反應的系統識別分析結果 ... 31 4.4.1 識別動態特性準確指標 ... 32 4.4.2 無噪訊與含噪訊數值模擬反應之識別結果 ... 32 4.4.3 濾波後對於時變系統的數值模擬反應識別結果之影響與比較 ... 33 第五章 結論與建議 ... 35 5.1 結論 ... 35 5.1.1 非時變系統 ... 35 5.1.2 時變系統 ... 35 5.2 建議 ... 36 參考文獻 ... 38 附表………. ... 43 附圖……. ... 55
viii
表目錄
表 4.1 七層樓剪力構架之非時變系統反應其理論值的系統識別結果 ... 43 表 4.2 七層樓剪力構架含噪訊之非時變系統反應經卡氏濾波過後之噪訊比 ... 44 表 4.3 七層樓剪力構架含 10%噪訊之訊號與經卡氏濾波後之系統識別結果... 45 表 4.4 七層樓剪力構架含 20%噪訊之訊號與經卡氏濾波後之系統識別結果... 46 表 4.5 七層樓剪力構架含 30%噪訊之訊號與經卡氏濾波後之系統識別結果... 47 表 4.6 七層樓剪力構架含 40%噪訊之訊號與經卡氏濾波後之系統識別結果... 48 表 4.7 七層樓剪力構架含噪訊之非時變系統反應經中值濾波過後之噪訊比 ... 49 表 4.8 七層樓剪力構架含 10%噪訊之訊號與經中值濾波後之系統識別結果... 49 表 4.9 七層樓剪力構架含 20%噪訊之訊號與經中值濾波後之系統識別結果... 50 表 4.10 七層樓剪力構架含 30%噪訊之訊號與經中值濾波後之系統識別結果... 50 表 4.11 七層樓剪力構架含 40%噪訊之訊號與中值濾波後之系統識別結果 ... 51 表 4.12 七層樓剪力構架含噪訊之非時變系統反應經離散小波濾波過後之噪訊比 ... 51 表 4.13 七層樓剪力構架含 10%噪訊之訊號與離散小波濾波後之系統識別結果... 52 表 4.14 七層樓剪力構架含 20%噪訊之訊號與離散小波濾波後之系統識別結果... 52 表 4.15 七層樓剪力構架含 30%噪訊之訊號與離散小波濾波後之系統識別結果... 53 表 4.16 七層樓剪力構架含 40%噪訊之訊號與離散小波濾波後之系統識別結果... 53 表 4.17 七層樓剪力構架含噪訊之非時變系統反應訊號經三種濾波方法處理後的系 識別結果綜合比較 ... 54ix
圖目錄
圖 2.1 含噪訊號的濾波流程 ... 55 圖 2.2 卡氏濾波器之模型示意圖 ... 55 圖 2.3 卡爾曼濾波器之分析流程 ... 56 圖 2.4 小波分析樹狀圖 ... 57 圖 2.5 小波重構示意圖 ... 58 圖 4.1 濾波與系統識別分析流程 ... 59 圖 4.2 a 七層樓之剪力構架示意圖 ... 60 圖 4.2 b 七層樓之剪力構架示意圖 ... 61 圖 4.3 a 七層樓剪力構架第一層樓之地震歷時反應圖(非時變系統) ... 62 圖 4.3 b 七層樓剪力構架第二層樓之地震歷時反應圖(非時變系統) ... 63 圖 4.3 c 七層樓剪力構架第三層樓之地震歷時反應圖(非時變系統) ... 64 圖 4.3 d 七層樓剪力構架第四層樓之地震歷時反應圖(非時變系統) ... 65 圖 4.3 e 七層樓剪力構架第五層樓之地震歷時反應圖(非時變系統) ... 66 圖 4.3 f 七層樓剪力構架第六層樓之地震歷時反應圖(非時變系統) ... 67 圖 4.3 g 七層樓剪力構架第七層樓之地震歷時反應(非時變系統) ... 68 圖 4.4 a 五層樓之剪力構架示意圖 ... 69 圖 4.4 b 五層樓之剪力構架示意圖 ... 70 圖 4.5 a 五層樓剪力構架第一層樓之地震歷時反應圖(時變系統) ... 71 圖 4.5 b 五層樓剪力構架第二層樓之地震歷時反應圖(時變系統) ... 72 圖 4.5 c 五層樓剪力構架第三層樓之地震歷時反應圖(時變系統) ... 73 圖 4.5 d 五層樓剪力構架第四層樓之地震歷時反應圖(時變系統) ... 74 圖 4.5 e 五層樓剪力構架第五層樓之地震歷時反應圖(時變系統) ... 75 圖 4.6 五層樓剪力構架無噪訊的時變系統反應之識別結果(理論值) ... 76 圖 4.7 五層樓剪力構架含 1%噪訊的時變系統反應之識別結果... 77 圖 4.8 五層樓剪力構架含 1%噪訊的時變系統反應經卡氏濾波器處理後之識別結果.... 78 圖 4.9 五層樓剪力構架含 1%噪訊的時變系統反應經中值濾波器處理後之識別結果.... 79 圖 4.10 五層樓剪力構架含 1%噪訊的時變系統反應經離散小波濾波處理後之識別結果 ... 80 圖 4.11 五層樓剪力構架含 3%噪訊的時變系統反應之識別結果 ... 81 圖 4.12 五層樓剪力構架含 3%噪訊的時變系統反應經卡氏濾波器處理後之識別結果.. 82x 圖 4.13 五層樓剪力構架含 3%噪訊的時變系統反應經中值濾波器處理後之識別結果.. 83 圖 4.14 五層樓剪力構架含 3%噪訊的時變系統反應經離散小波濾波處理後之識別結果 ... 84 圖 4.15 五層樓剪力構架含 5%噪訊的時變系統反應之識別結果... 85 圖 4.16 五層樓剪力構架含 5%噪訊的時變系統反應經卡氏濾波器處理後之識別結果.. 86 圖 4.17 五層樓剪力構架含 5%噪訊的時變系統反應經中值濾波器處理後之識別結果.. 87 圖 4.18 五層樓剪力構架含 5%噪訊的時變系統反應經離散小波濾波處理後之識別結果 ... 88
1
第一章 緒論
台灣位處歐亞板塊與菲律賓板塊之環太平洋之地震帶，地震相當頻繁，尤其在民 國八十八年「九二一大地震」發生之後，許多民眾不僅在經濟上受到重大損失，甚至 失去寶貴性命，造成許多社會問題，也因此使得國人日益重視結構物耐震的設計與以 及安全性。 土木工程系統最大的問題在於結構物實際上表現出來的力學行為與理論上的力學 行為差異較大且相當複雜，加上現地施工品質不易掌握，因此藉由現地量測的資料來 識別結構物的動態特性以及其健康狀況受到相當程度地重視。 「九二一大地震」發生之後，加速隔減震裝置越來越廣泛地運用於補強或修復結 構物，而許多近期研發之「半主動控制元件」例如：半主動摩擦組尼器、半主動控制 磁流變消能器、蓄壓式油壓半主動減震器以及電流變阻尼器……等等，亦相繼地被納 入於新開發結構物之設計中。當此類半主動控制元件被裝置於結構物上，該結構系統 受到外力作用時，其勁度和阻尼將會隨時間而改變，此即為系統識別領域中所稱之時 變系統( Time-varient system)。現今科技日新月異，吾人可以透過系統識別( System Identification)技術來處理相關的工程問題，系統識別技術可以識別結構系統之動態特 性並協助建立非線性模型，並幫助檢驗半主動控制元件等非線性元件的應用是否完善， 評估結構系統的健康狀況。 另外，若是結構系統在受外力作用之時表現出材料的非線性行為(Non- linear)，亦 可視其為非線性時變系統，吾人可透過系統識別評估該結構系統隨時間改變之動態特 性，瞭解其力學與其動態行為之影響。1.1 前言
在工程應用上，大部分碰到的問題為時變系統，可用來描述非穩態過程中之各種 特性，TVARX(Time Varying Autoregressive with exogenous input)模型常被用來表示時 變系統之輸入與輸出之關係。如前述所知，發展出一準確性高且具有高效率的系統識別技術，對於土木工程領 域的應用具有相當大地重要性，然而不論是時變系統還是非時變系統，普遍存在「雜 訊」(或稱噪訊)幹擾的問題，在系統識別過程中，噪訊幹擾會嚴重影響到系統識別的 結果。因此，利用「濾波」(filtering，又稱去噪，de-noising)方法，盡可能地降低噪訊
2 對系統識別的結果之影響，成為在系統識別領域中一個相當重要的課題。 目前在訊號處理領域中，所發展出來的濾波方法有非常多種；針對不同類型的訊 號，亦發展出各種不同的濾波方法，每一種濾波方法都有其優缺點，也都有其適合的 訊號類型。濾掉太多噪訊易導致原始的訊號失真，喪失原有的訊號特性；濾掉太少則 殘留的噪訊太多，使得系統識別的結果大受影響。 本研究針對時變系統以及非時變線性系統訊號，利用目前使用上較為普遍的三種 濾波方法，分別是卡爾曼濾波演算法、中位值濾波法以及離散小波濾波法，進行其對 系統識別結果影響之比較與討論。
1.2 文獻回顧
(1) 卡爾曼濾波演算法( Kalman filtering algorithm，亦稱卡氏濾波器)：1960 年起， Kalman 發表了一種濾波的演算法[1,2]，該濾波方法並以他命名，此後卡氏濾波器 在濾波領域逐漸受到重視並持續發展，傳統的卡氏濾波器為基本型卡爾曼濾波器 (The Basic Kalman filter)，也是本研究中所使用的濾波方法。多數研究指出，該方 法在處理離散系統線性濾波的問題中可以得到良好的效果，目前主要應用於自主 或輔助定位導航、測量、信號處理、以及自動控制等領域；此外，Gelb(1974) [3]、 Lewis(1986) [4]、Brown(1992)[5]、Jacobs(1993)[6]等人對卡氏濾波器做了更深入地 研究；Harvey (1989)[7]、Welch and Bishop(2001)[8]等人也對其做了詳盡的介紹。 至於用於非線性問題時，一般則使用擴展型卡爾曼濾波器(The Extended Kalman Filter，簡稱 EKF)：Ljung(1979)[9]；Hoshiya and Saito(1984)[10]；Roman(1987)[11]； Julier(1997)[12]；Vander Merwe and Wan(2001)[13]；張(2008)[14]。
(2) 中位值濾波法( Median filter，亦稱中值濾波器)：1983 年，Bednar 提出時間方向的 中值濾波，將中值濾波法引入地震探勘領域[15]；1984 年，Brownrigg 發表了考慮 權重參數的中值濾波器[16]；1995 年，Duncan 也陸續發表了不少關於中值濾波應 用於地震探勘處理中的論文[17,18]；2006 年，Liu 等利用二維多級中值濾波技術進 行白噪消除[19]，並於 2009 年發表了非穩態時變中值濾波技術[20]；中值濾波器 (median filter)及其相關濾波器在各領域都有得到廣泛的發展：Bednar(1984)[21]， Lee and Kassam(1985) [22]；Fitch, et al.(1986)[23]。
(3) 離散小波去噪(Discrete wavelet de-noising)：1992 年，Mallat 等人提出了基於信號 的奇異性，利用 Lipschitz 指數在多尺度上對信號、圖像和噪聲的數學特性進行描
3 述，稱模極大值重構濾波方法[24,25]；Lu 等人直接結合小波變換理論與傳統的多 尺度信號處理方法，提出了一種效果優良的小波濾波法[26,27]；Simoncelli, et al. 提出利用在小波域對噪聲的小波係數進行抑制達到去噪的效果[28]；1994 年， Moulin 在功率譜之估計問題的研究上得到與閾值演算法相類似的結果[29]；同年， Dohono 正式定義了「去噪」，通過對小波係數閾值化處理而去除噪聲，利用均方 差來做為最優去噪效果的指標，並提出了小波域閾值濾波演算法，並有大量的理 論和成果[30~35]；Kirm 等人利用最小描述長度( Minimum description length,MDL) 準則，得到了相同的閾值[36]；Moulin 在研究譜估計方面的問題提出了另一種選 擇閾值的方法[37]。
1.3 研究內容
本研究分為兩大部分進行，非時變線性系統以及線性時變系統訊號的濾波處理和 系統識別。由於噪訊對系統識別之結果影響很大，然雜訊在實際量測反應上卻又是不 可避免的，因此本研究在非時變系統部分，利用六層樓及七層樓構架之數值模擬地震 反應，分別加上 10%、20%、30%、40%的高斯白噪聲(Gaussian white noise，簡稱白噪 聲或白噪)，經過前述三種濾波方法濾波之後，討論其濾波之後對系統識別結果的影響； 而在線性時變系統部分，則利用五層樓構架之數值模擬地震反應，分別加上 1%、3%、 5%的白噪，經過前述三種濾波方法濾波之後，探討其濾波之後對系統識別結果的比較 與影響。 論文架構如下： 第一章 緒論：主要內容為文獻回顧以及研究內容之簡介。 第二章 濾波方法理論與介紹：分別介紹本研究中所使用之三種濾波方法的基本理論與 操作方法。 第三章 架構 TVARX 模型以及 ARX 模型：介紹如何建立該兩種模型與其基礎理論。 第四章 數值模擬之系統識別結果的分析：分為非時變系統和時變系統兩大部分來討論 其系統識別結果。 第五章 結論與建議：總結此研究中非時變系統與時變系統之濾波效果與系統識別分析 結果。4
第二章 濾波方法介紹
2.1 前言
噪訊在實際應用上難以避免，卻會對系統識別的結果造成一定程度的影響。為減 少噪訊的在系統識別應用中所造成的問題，本研究擬利用三種濾波方法進行訊號處理， 分別是卡爾曼濾波法、離散小波濾波法以及中位值濾波法(請參閱 1.2 節之文獻回顧)， 希望含噪訊號可以經由濾波過後，降低雜訊對於系統識別結果的影響。如圖 2.1 所示， 理論上我們希望可以藉由濾波器的處理將噪訊的含量及影響減至最低，並盡可能保留 真實的訊號(即盡量減低真實訊號在濾波過程中的損失)。上述這三種濾波法都是可以 用來處理線性及非線性系統的濾波方法，以下將對這三種濾波方法的原理進行基本的 介紹：2.2 節介紹卡爾曼濾波器的理論及應用；2.3 節介紹中位值濾波法(亦稱中值濾波 器)；2.4 節介紹離散小波濾波的理論及應用。2.2 卡爾曼濾波演算法之介紹
1960 年，卡爾曼(Rudolph E. Kalman)發表了以遞歸方法解決離散數據線性濾波問 題的論文，故此濾波方法以卡爾曼命名。卡爾曼濾波(Kalman filter；以下簡稱卡氏濾 波器) 由一系列的遞歸數學公式描述，是一種高效率的遞歸濾波器(亦稱自回歸濾波器) ，它能夠從一系列的不完全及包含噪訊的測量(measurement)中，估計動態系統的狀態。 一般而言，量測點的位移、速度、加速度的測量值往往都會有噪訊，卡氏濾波器利用 量測點的動態資訊，從一組包含噪聲的序列，設法除去噪聲的影響，預測其位移及速 度，良好地估計出量測點的狀態。 簡單的說，就是以最小均方誤差為估計的最佳準則，來尋求一套遞迴估算的演算 法，當輸入為高斯白噪聲(Gaussian white noise，簡稱白噪，即常態分佈的隨機噪訊)時， 使得期望的量測值和實際的輸出值之間的均方根誤差達到最小。採用訊號與噪訊的狀 態空間模型，利用前一時刻的估計值和當前時刻的觀測值來更新對狀態變數的估計， 求出當前時刻的估計值。卡氏濾波器在許多工程應用(如雷達、動態定位系統、GPS、經濟學、時間序列模 型……等)中都有良好的發展亦是控制理論以及控制系統工程中的一個重要課題。
5
2.2.1 卡氏濾波器的基本動態系統模型
卡氏濾波器建立在隱馬爾可夫模型(Hidden Markov Model)上，其模型假設當 t 時刻
的真實狀態是從(t − 1)時刻的狀態演化而來，可用下式表示之：
𝒙𝑡 = 𝑭𝑡𝒙𝑡−1+ 𝑩𝑡𝒖𝒕+ 𝒘𝑡 (2.1) 其中𝑭𝑡為作用在𝒙𝑡−1的狀態變換模型；𝑩𝑡為作用在控制器向量𝒖𝑡的輸入控制模型；𝒘𝑡 為符合均值為零的過程噪聲(process noise)且協方差矩陣為𝑸𝑡的多元常態分佈，𝒘𝑡~ N(0, 𝑸𝑡)。 對於一個真實狀態𝐱𝑡的測量𝒛𝑡，在時刻為 t 時滿足下式： 𝒛𝑡= 𝑯𝑡𝒙𝑡+ 𝒗𝑡 (2.2) 其中𝑯𝑡為觀測模型，本研究在濾波演算過程中將其設為𝑯𝑡= [ 1 0 ]；而 𝒗𝑡 為均值為 零的觀測噪聲(observation noise)，且協方差矩陣為𝑹𝑡並滿足𝒗𝑡~N(0, 𝑹𝑡)。假設𝒘𝑡與𝒗𝑡 皆為互相獨立且正態分佈的白色噪聲；實際的應用中， 𝒘𝑡與𝒗𝑡通常為有色噪聲，兩者 亦可能存在互相關的情況，而𝑸𝑡與𝑹𝑡則可能會隨著疊代計算而產生變化，但在這裡我 們將𝑸𝑡與𝑹𝑡假設為常數[8]。在本研究中，我們假設𝑸𝑡為 0.001；且在此稱 GR 值為「猜 測的噪訊值」，在濾波之前，先假設 GR 值分別為 10%、20%、30%、40%，其中 GR 值即代表式中的 ，濾波過程中我們將 假設為 0.1、0.2、0.3、0.4，欲觀察猜測的噪 訊含量是否與實際上的噪訊含量一致，而此過程對於過濾噪訊與系統識別的結果又將 影響與否。 卡氏濾波器模型請參考圖 2.2，其中圓圈代表向量，方塊代表矩陣，星號代表高斯 噪聲，其協方差矩陣在右下方標示。
2.2.2 卡氏濾波器的推導及運算
在線性的條件之下，卡氏濾波器的狀態可由以下兩個變數表示： 𝒙̂ 𝑡 | 𝑡：當時刻為 t 時的狀態估算， 𝑷 𝑡 | 𝑡：誤差的協方差矩陣(covariance matrix)，用來度量估計值的精確程度。
6 基本型的卡氏濾波器的運算包含預測及更新兩個階段： (1) 預測：依式(2.1)，利用前一刻狀態的估算，做出對當前時刻狀態的估算 預測狀態： 𝒙̂ 𝑡 | 𝑡−1 = 𝑭𝑡𝒙̂ 𝑡−1 | 𝑡−1+ 𝑩𝑡𝒖𝒕−𝟏 (2.3) 預測估計協方差矩陣： 𝑷 𝑡 | 𝑡−1 = 𝑭𝑡𝑷 𝑡−1 | 𝑡−1𝑭𝑡𝑇+ 𝑸𝑡 (2.4) (2) 更新： 𝒙̂ 𝑡 |𝑡 = 𝒙̂ 𝑡 | 𝑡−1+ 𝑲𝑡𝒚̃𝑡 (2.5) 𝑷𝑡|𝑡= (𝑰 − 𝑲𝑡𝑯𝑡)𝑷𝑡|𝑡−1 (2.6) 其中， 𝒚̃𝑡 = 𝒛𝑡− 𝑯𝑡𝒙̂𝑡|𝑡−1 (2.7) 𝑲𝑡 = 𝑷𝑡|𝑡−1𝑯𝑡𝑇𝑺𝑡−1 (2.8) = | −1 + (2.9) 𝒚̃𝑡為測量餘量(measurement residual)；𝑺𝑡為測量餘量的協方差矩陣；𝑲𝑡為卡爾曼增益； 𝒙̂ 𝑡 |𝑡為更新的狀態估算；𝑷𝑡|𝑡為更新的協方差矩陣之估算，只適用於最優卡爾曼增益。 假設模型準確，且𝒙̂0|0和𝑷0|0反映了最初狀態的分佈，那麼可知協方差矩陣亦準確 地反映了估算的協方差，即 𝑷𝑡|𝑡= 𝑐𝑜𝑣(𝒙𝑡− 𝒙̂ 𝑡|𝑡) (2.10) 𝑷𝑡|𝑡−1= 𝑐𝑜𝑣(𝒙𝑡− 𝒙̂ 𝑡|𝑡−1) (2.11) 𝑺𝑡= 𝑐𝑜𝑣(𝒚̃𝑡) (2.12) 其中𝑐𝑜𝑣(𝑎) = [𝑎𝑎𝑇]， [ ]代表 的期望值。 欲推導後驗( posteriori)協方差矩陣，首先將式(2.5)代入式(2.10)，可得： 𝑷𝑡|𝑡= 𝑐𝑜𝑣[𝒙𝑡− (𝒙̂ 𝑡 | 𝑡−1+ 𝑲𝑡𝒚̃𝑡)] (2.13)
7 再將式(2.2)及(2.7)代入式(2.13)，得： 𝑷𝑡|𝑡= 𝑐𝑜𝑣{𝒙𝑡− [𝒙̂ 𝑡 | 𝑡−1+ 𝑲𝑡(𝑯𝑡𝒙𝑡+ 𝒗𝑡− 𝑯𝑡𝒙̂𝑡|𝑡−1)]} (2.14) 整理可得： 𝑷𝑡|𝑡= 𝑐𝑜𝑣[(𝑰 − 𝑲𝑡𝑯𝑡)(𝒙𝑡− 𝒙̂𝑡|𝑡−1) − 𝑲𝑡𝒗𝑡] (2.15) 由於測量誤差𝒗𝑡與其他項不相關，因此(2.13)式可改寫為： 𝑷𝑡|𝑡= 𝑐𝑜𝑣[(𝑰 − 𝑲𝑡𝑯𝑡)(𝒙𝑡− 𝒙̂𝑡|𝑡−1)] + 𝑐𝑜𝑣(𝑲𝑡𝒗𝑡) (2.16) 利用協方差矩陣的性質，式(2.14)可以改寫為： 𝑷𝑡|𝑡= (𝑰 − 𝑲𝑡𝑯𝑡)𝑐𝑜𝑣(𝒙𝑡− 𝒙̂𝑡|𝑡−1)(𝑰 − 𝑲𝑡𝑯𝑡)𝑇+ 𝑲 𝑡𝑐𝑜𝑣(𝒗𝑡)𝑲𝑡𝑇 (2.17) 最後，將式(2.11)代入(2.15)，並令𝑹𝑡 = 𝑐𝑜𝑣(𝒗𝑡)，整理可得下式： 𝑷𝑡|𝑡= (𝑰 − 𝑲𝑡𝑯𝑡)𝑷𝑡|𝑡−1(𝑰 − 𝑲𝑡𝑯𝑡)𝑇+ 𝑲𝑡𝑹𝑡𝑲𝑡𝑇 (2.18) 上式對於任一卡爾曼增益𝑲𝑡皆成立。卡氏濾波器的分析流程圖請參考圖 2.3。
2.3 中位值濾波法之介紹
中位值濾波法( median filter；簡稱中值濾波器)是基於排序統計理論中的一種 可以有效地抑制噪訊的非線性信號處理技術。中值濾波器能夠有效地消除非平穩 信號中的峰值噪聲， 經常被使用於去除圖像的噪訊(避免邊緣模糊或進行邊緣檢測) 或者其它信號中的噪訊，對於消除斑點噪聲(speckle noise，又稱散斑幹擾或散斑噪聲， 為一種顆粒狀的雜訊，通常由反射信號的隨機波動所造成，會增加二維圖像局部區域 的平均灰階，造成圖像判讀困難)和椒鹽噪聲(salt-and-pepper noise，為隨機產生的黑點 和白點狀，通常由於訊號的快速瞬變或故障切換時產生)效果尤其顯著[38]。 中值濾波器的基本原理是把信號序列中某一點的值，用該點附近的一個鄰域 中各個點值由大排列到小，取其中的中位數來替代該點本來的值，讓周圍的圖元 值接近的真實值，從而消除孤立的噪聲點。目前，中值濾波器運用 在地震資料處 理方面亦有相當不錯的成果[39~42]，本研究將中值濾波器應用於一維數值模擬8 多自由度的剪力構架受地震反應之訊號濾波與系統識別，其結果將在第四章討 論。
2.3.1 中位值濾波法之公式推導
中值濾波器是一種能有效抑制噪訊的非線性濾波技術，基本原理是將一個序 列中的某值用涵蓋該點的跨度(N)中各個點值的中位數值取代，藉此消除孤立的 噪聲點，使得誤差的絕對值之和達到最小值來估算參數的方法。其推導過程如 下： 令一組觀測值{𝑥𝑖}，其中 1≤ i ≤ N，其絕對誤差和表示為： Q = ∑|𝑥̅ − 𝑥𝑖| 𝑁 𝑖=1 (2.19) 欲使 Q 達到最小，必須令dQ / d𝑥̅ = 0 dQ d𝑥̅ = 𝑑 𝑑𝑥̅(∑|𝑥̅ − 𝑥𝑖| 𝑁 𝑖=1 ) = 𝑑 𝑑𝑥{∑[(𝑥̅ − 𝑥𝑖)2] 1 2 𝑁 𝑖=1 } = ∑ 𝑥̅ − 𝑥𝑖 |𝑥̅ − 𝑥𝑖| 𝑁 𝑖=1 = ∑ 𝑠𝑖𝑔𝑛(𝑥̅ − 𝑥𝑖) 𝑁 𝑖=1 = 0 (2.20) 其中sign(ω𝑗,𝑖)為符號函數： sign(ω𝑗,𝑖) = { −1，ω𝑗,𝑖 < 0 0， ω𝑗,𝑖= 0 1，ω𝑗,𝑖 > 0 (2.20a) 由式(2.20)可知，若要滿足 ∑𝑁𝑖=1𝑠𝑖𝑔𝑛(𝑥̅ − 𝑥𝑖) = 0，所選取之 𝑥̅ 必須使得 N 個 𝑥𝑖 中有 N/2 個 𝑥𝑖大於 𝑥̅，即 𝑥̅ 正好為{𝑥𝑖}序列中的值按大至小順序排列時的中間 值，此即稱為中位數；一般而言，取 N 為奇數，若該序列的個數 N 為偶數，則 取最中間兩個數的平均值為該序列的中位數(亦稱為中值)。 中值濾波器的操作過程分為兩個步驟： (1) 從訊號中的某個採樣視窗取出 N 個數據由小到大進行排序。9 (2) 取該序列的中位數取代待處理的數據即可。 此外，中值濾波器的效果取決於「跨度」(N)：每個訊號濾波所適用的跨度並 不一定，跨度太小，信號不夠平滑，達不到去噪效果；跨度太大，平滑效應會導 致有效信號會被壓制。因此，如何選擇適當的跨度是中值濾波的關鍵，本研究所 使用之跨度 N 約為 200~250 之間。
2.4 離散小波轉換之介紹
小波變換的概念是由法國工程師 Morlet 在 1974 年首先提出的，它是一個時間和 頻率的局域變換，能夠有效地從訊號中提取資訊，透過伸縮或平移等運算功能，構成 𝐿2(𝑅)函數空間的一組完整的基底函數，對函數或信號進行多尺度分析。小波分析為一 種強大的數學工具，目前被廣泛地應用於信號處理的領域，可用一簡單的多分辨分析 小波樹狀圖來表示(圖 2.4)。1988 年 Mallet 提出了多分辨分析(Multi-Resolution Analysis，MRA)的概念，將所 有的正交小波基的構造法統一起來，說明瞭小波的多解析度特性。小波分析不同於短 時傅立葉變換在單解析度上的缺點，具有多分辨分析的優點，在時域和頻域都可以處 理訊號的局部訊息。小波變換具有下列特性：時頻局部化特性、多解析度特性以及小 波基選擇的靈活性。
2.4.1 離散小波變換的推導
離散小波變換( Discrete Wavelet Transform)在數值分析和時頻分析中很有用。離散 小波轉換顧名思義就是離散的輸入以及離散的輸出。 在連續小波中，考慮小波函數： 𝛹𝑎,𝑏(𝑡) = |𝑎| − 1 2 𝛹(𝑡−𝑏 𝑎 ) (2.21) 其中𝑏 ∈ 𝑅，𝑎 ∈ 𝑅+，且𝑎 ≠ 0。 在實際運用上，由於數據處理上的問題，連續小波必須經過離散化。限制 a >0， 所得之相容性條件變為：
10 𝐶𝛹 = ∫ |𝛹̂ (𝜔̅ )| 2 (𝜔̅ ) ∞ 0 𝑑𝜔 ̅̅̅ < ∞ (2.22) Meyer 於 1986 年發表了具有一定衰減特性之光滑小波函數，此函數之二進制尺度 與平移足以表示𝐿2(𝑅)空間中之所有函數。將連續小波變換中的尺度因數 a 和平移因數 b 進行二進位之離散化，取 𝑎 = 𝑎0𝑗，𝑏 = 𝑘𝑎0𝑗𝑏0，其中 𝑘 , 𝑗 ∈ 𝑍。 可得離散小波函數為：
𝛹𝑗,𝑘(𝑡) = 𝑎0 − 𝑗 2 𝛹(𝑡−𝑘𝑎0 𝑗𝑏 0 𝑎0𝑗 ) = 𝑎0 − 𝑗2 𝛹(𝑎 0−𝑗− 𝑘𝑏0) (2.23) 通常我們對其進行二進制離散化，取𝑎0 =2 且 𝑏0 =1，則𝛹𝑗,𝑘(𝑡)可表示為： 𝛹𝑗,𝑘(𝑡) = 2 − 𝑗 2 𝛹(2−𝑚𝑡 − 𝑛) (2.24) 此亦稱為二進制小波(Dyadic Wavelet)。 我們將離散化小波變換表示為： 𝐷𝑊𝑓(𝑗, 𝑘) = ∫ 𝑓(𝑡)𝛹𝑅 ̅𝑗,𝑘(𝑡)𝑑𝑡 (2.25) 若用內積的方式表示則寫為： 𝑫𝑾𝑓(𝑗, 𝑘) = 〈 𝒇(𝑡), 𝜳𝑗,𝑘(𝑡) 〉 𝑗, 𝑘 ∈ 𝑍 (2.26) 值得注意的是，離散小波轉換𝐷𝑊𝑗,𝑘基本上仍然將信號𝑓(𝑡)經過一系列帶通濾波器而輸 出，不同的是其帶通濾波器之中心頻率和帶寬由於 a 之離散取樣而成為一系列之離散 值，並且濾波之後的輸出也因為 b 的離散取樣而成為若干離散取樣值。 若取具正交性之小波： 〈𝛹𝑚,𝑛，𝛹𝑗,𝑘〉 = {1, 當𝑚 = 𝑗 且 𝑛 = 𝑘 0, 當𝑚 ≠ 𝑗 或 𝑛 ≠ 𝑘 ， (2.27)
11 則對任一函數𝑓(𝑡) ∈ 𝐿2(𝑅)，可經由小波轉換表示成： 𝑓(𝑡) = ∑𝑁𝑛=1𝐶𝑉(𝑀, 𝑛)Φ𝑀,𝑛+ ∑𝑀𝑚=1∑𝑁𝑛=1𝐶𝑊(𝑚, 𝑛)𝛹𝑚,𝑛 (2.28) 其中Φ代表尺度函數(Scaling function)。 當取不同之 m 時，即是將𝐿2(𝑅)分解成不同的空間，各個不同分解空間之關係可 以為(Mallat,1989a and 1989b)： 𝐿2(𝑅) = 𝑊1⨁ 𝑉1 = 𝑊1⨁ 𝑊2⨁ 𝑉2 = 𝑊1⨁ 𝑊2⨁ ⋯ ⨁𝑊𝑀⨁ 𝑉𝑀 (2.29) 其中 𝑊𝑖⨁ 𝑉𝑖 且 𝑉𝑖 = 𝑉𝑖+1⨁𝑊𝑖+1，而 𝑉𝑖空間之訊號較𝑊𝑖低頻。𝛹𝑚,𝑛與Φ𝑚,𝑛分別代表𝑊𝑚 與 𝑉𝑚之正交基底。理論上，所選取之 M 足夠大時， 𝑉𝑀將趨近於零空間。
2.4.2 離散小波在濾波上的應用
信號和噪訊在小波域中呈現不同的狀態，兩者的小波係數會隨著尺度的不同而產 生不同的變化，隨著尺度分解的階數越多，噪聲係數的幅值會很快地衰減，而原始訊 號係數的幅值基本上維持不變。濾波的重點就在於盡可能減小由噪訊所產生的係數， 且盡可能地保留原始訊號的係數，最後將經過處理的小波係數重構回去，得到原始信 號估計的最優化。 小波濾波的三大基本假設： (1) 噪訊經小波變換之後，大部分的小波係數趨近零。 (2) 噪訊均勻地分佈在所有係數之中。 (3) 噪訊比不至於過高[30]。 小波的信號濾波( filtering，亦可稱為去噪，de-noising)的過程可概略分為三個步 驟： (1) 選擇一個小波分解的尺度 (階數)及小波基，進行分解計算。 (2) 選擇一個閾值(threshold)，對各個分解尺度下的小波係數進行非線性處理。 (3) 將處理過後的低頻及高頻係數進行重構(小波逆變換)，見圖 2.5。 普遍來說，地震訊號中的噪訊只存在於前三或四尺度的分解結果中[43]，但並無特定 的判斷準則來確定小波分解的尺度，只能大概估算。12 小波基的選擇又是另一門學問，一個良好的濾波效果，需要在緊支撐性、平滑性 以及對稱性等達到一定的水準。然而目前階段而言，這三點不太可能同時達到滿足： 緊支撐性與平滑性不可兼得，要有較好的平滑特性，則緊支撐長度需要增加；反之， 緊支撐的長度較小，雖然可以保留訊號的局部特性，但是相對的平滑性也會下降。本 文中，我們參照柳等人[43]的建議來選取最優小波基，採用 sym4 小波基。 在小波濾波的過程中，閾值對於濾波的效果具有決定性的影響，目前有許多文獻 提出了各種閾值的選取方法及理論，例如：通用閾值法是 Dohono 和 Johnstone 將多維 獨立正態變量決策理論應用於高斯噪聲模型所得出的通用閾值t = σ√2lnN[32]；潘等 人驗證了基於正態分佈的自適應閾值 𝑡 = 𝑐𝜎，𝑐 ∈ [3.0,4.0][44,45]。極小化風險閾值法
則利用SURE(Stein’s Unbiased Risk Estimator)法[33,46]、交叉驗證(Cross-Validation, CV)
演算法[47]、廣義交叉驗證(Generalized Cross- Validation, GCV)演算法[48,49]對均方差 函數進行估計以得到濾波後的結果與其真實訊號之間的差異來確定閾值……等等。不 過，每個訊號的去噪效果對閾值的反應不一，目前並沒有很明確的方法可以判別何種 閾值會產生比較好的去噪效果。一般情況下，閾值的選定會採取試誤法(trial and error method)。 閾值確定後，再對小波係數進行閾值量化處理，目前主要的方法有硬閾值法及軟 閾值法，其基本概念都是去除小幅值的係數，對較大幅值的係數進行收縮或保留 [30~35]。硬閾值法因為有不連續性，通常會使得結果產生較大的方差，優點為可以保 留信號邊緣的局部特徵，缺點為易造成偽吉布斯(Gibbs)現象。軟閾值法因為收縮性質 較大，濾波結果相對平滑許多，卻容易造成濾波結果有較大的偏差，且使得邊緣模糊 等失真現象。 硬閾值法的基本處理方式為當某時刻的小波變換值大於閾值 t 時將保留原值，否 則置零，其公式為： 𝜃̂𝑗,𝑖 = {0 ， |ω𝑗,𝑖| ≤ 𝑡 ω𝑗,𝑖 ，|ω𝑗,𝑖| ≤ 𝑡 ， 其中 𝑡 = σ√2𝑙𝑛𝑁 (2.30) 軟閾值法則是當某時刻的小波變換值大於閾值 t 時，將其值減去 t 值，反之則將 其設為零，公式如下： 𝜃̂𝑗,𝑖 = sign(ω𝑗,𝑖) ∙ ( |ω𝑗,𝑖| − 𝑡) ， 𝑡 = σ√2𝑙𝑛𝑁 (2.31) 其中sign(ω𝑗,𝑖)為符號函數，定義同式(2.20a)。
13
1995 年，Gao 提出了結合硬閾值法與軟閾值法的一種折衷形式，稱為半軟閾值法 (semisoft shrinkage) [50]，1997 年並推導出基於該方法的 Minimax 閾值[51]，至於上閾 值和下閾值的選取方式，請參閱文獻[50]，於此不再贅述。在此法中，需要確定兩個 閾值，會增加計算上的複雜度，但可以截長補短，保留硬閾值法及軟閾值法各自的優 點。半軟閾值法的表示式如下： 𝜃̂𝑗,𝑖 = { 0 ，|ω𝑗,𝑖| ≤ 𝑡2 sign(ω𝑗,𝑖)𝑡2( |ω𝑡 𝑗,𝑖|−𝑡1) 2−𝑡1 ， 𝑡1 < |ω𝑗,𝑖| ≤ 𝑡2 ω𝑗,𝑖 ，|ω𝑗,𝑖| ≤ 𝑡2 (2.32) 其中𝑡1為上閾值，𝑡2為下閾值。
14
第三章 架構 ARX 模型以及 TVARX 模型
3.1 前言
本章將介紹兩種系統識別模型。第一個是介紹如何透過連續小波轉換，將量測之 訊號轉換至小波域中，並於小波域中架構量測系統之 ARX 模型(Autoregressive with exogenous input Model)，再進一步估算系統模態。第二個是移動最小平方差法( Moving Least Square method，MLS)架構 TVARX 模型( Time Varying Autoregressive with exogenous input Model)，識別系統之瞬時模態特性。此方法主要是透過數據之輸入與 輸出反應架構 TVARX 模型，再進一步估計系統之瞬時模態參數。 3.2 節介紹連續小波變換之理論基礎；3.3 節介紹如何於小波域中架構 ARX 模型， 並估算結構系統動態特性；3.4 節介紹以移動最小平方差法架構 TVARX 模型；3.5 節 介紹瞬時模態估算的方法。
3.2 連續小波變換之理論基礎
小波變換為近十幾年來在信號處理領域中一個熱門的研究主題，目前廣泛地被運 用在圖像壓縮、特徵提取、信號濾波等方面。小波變換是傅立葉(Fourier)變換的延伸， 傳統的傅立葉變換因為無法描述信號的局部頻率特徵，因此不再是分析非平穩信號的 有效工具。小波變換既保留了傅立葉變換的優點，亦彌補了傅立葉變換在信號分析方 面的不足；因此，小波變換可視為傅立葉變換之改良。 給定一函數 𝑓(𝑡) ∈ 𝐿2(𝑅)，其連續小波轉換定義如下[52]： 𝑊𝑓(𝑎, 𝑏) = 1 √𝑎∫ 𝑓(𝑡) 𝛹̅( 𝑡−𝑏 𝑎 ) 𝑅 𝑑𝑡 = 〈 𝑓, 𝛹𝑎,𝑏 〉 (3.1) 其中𝛹(𝑡)為母小波函數(mother wavelet)，𝛹̅(𝑡)為𝛹(𝑡)的共軛函數，𝛹𝑎,𝑏為小波函數， 乃是經由母小波函數𝛹(𝑡)伸縮和平移所形成之基底函數，其定義為： 𝛹𝑎,𝑏(𝑡) = |𝑎| − 1 2 𝛹(𝑡−𝑏 𝑎 ) (3.2)15 小波逆變換的定義為： 𝑓(𝑡) = 𝐶1 𝛹∫ ∫ 1 𝑎2 ∞ −∞ ∞ −∞ 𝑊𝑓(𝑎, 𝑏)𝛹( 𝑡−𝑏 𝑎 )𝑑𝑎𝑑𝑏 (3.3) 𝐶𝛹的定義為： 𝐶𝛹 = ∫ |𝛹̂ (𝜔)| 2 (𝜔) 𝑅 𝑑𝜔 < ∞ (3.3a) 其中𝛹̂(𝜔)為𝛹(𝑡)之傅立葉變換。 有別於傅立葉變換中之基底函數被固定為𝑒𝑖𝜔𝑡，小波變換之母小波函數則有多種 選擇性，但要滿足以下兩個條件： (1) 小波函數必須為平方可積(Square Integrable)或函數具有有限能量。 (2) 小波函數及其各階倒數在 t → ∞ 處急遽下降收斂至零，即為有限支撐長度。 依式(3.1)定義，若母小波函數
t 在時間域，可將之視為中心為t，半徑為之 窗函數，𝛹𝑎,𝑏之中心在b + 𝑎t∗，且半徑為𝑎Δ。式(3.1)表示之小波轉換，只取時間窗： [𝑏 + 𝑎𝑡∗− 𝑎Δ , 𝑏 + 𝑎𝑡∗+ 𝑎Δ ] (3.4) 內𝑓(𝑡)的局部訊息。式(3.4)之窗中心t及半徑之定義分別為： 𝑡∗= 1 ‖𝑓‖2∫ 𝑡|𝑓(𝑡)|2 ∞ −∞ 𝑑𝑡 (3.5) ∆= 1 ‖𝑓‖{∫ (𝑡 − 𝑡∗)|𝑓(𝑡)|2 ∞ −∞ 𝑑𝑡} 1 2 (3.6) 其中‖𝑓‖2 = ∫ |𝑓(𝑡)|∞ 2 −∞ 𝑑𝑡。 由於 1 2𝜋𝛹𝑎,𝑏(𝜔) = 𝑎𝑒𝑖𝑏𝜔 2𝜋√|𝑎| 𝛹̂(𝑎𝜔) (3.7)16 對式(3.1)進行傅立葉變換，可得： |𝑊̂𝑓(𝑎, 𝜔)| = √𝑎|𝛹̂(𝑎𝜔)||𝑓̂(𝜔)| (3.8) 其中𝑊̂𝑓(𝑎, 𝜔)代表𝑊𝑓(𝑎, 𝑏)之傅立葉變換，根據上式可定義𝑓之頻率窗範圍為[53]： [𝜔 ∗ 𝑎 − ∆̂ 𝑎 , 𝜔∗ 𝑎 + ∆̂ 𝑎] (3.9) 其中 𝜔∗ 與 ∆̂ 為在頻率域計算𝛹̂之中心與半徑，定義如下： 𝜔∗ = 1 ‖𝑓̂‖2∫ 𝜔|𝑓̂(𝜔)| 2 ∞ −∞ 𝑑𝜔 (3.10) ∆̂= ‖𝑓̂‖1 {∫ (𝜔 − 𝜔−∞∞ ∗)|𝑓̂(𝜔)|2𝑑𝜔} 1 2 (3.11) 式(3.9)表示之窗具有中心頻率 𝜔𝑎∗ 且帶寬為 2∆̂𝑎 ，利用式(3.24)之小波變換，其將提 供一時間頻率窗： [𝑏 + 𝑎𝑡∗− 𝑎Δ , 𝑏 + 𝑎𝑡∗+ 𝑎Δ ] × [𝜔∗ 𝑎 − ∆̂ 𝑎 , 𝜔∗ 𝑎 + ∆̂ 𝑎] (3.12) 𝛹𝑎,𝑏(𝑡)所定義之時間－頻率窗口之矩形形狀，隨伸縮因數 a 變動而改變。當 a 值 大時，時間視窗較寬而頻率視窗變窄，適合分析低頻信號；而當 a 值小時，時間視窗 較窄而頻率視窗寬，適合偵測高頻現象，所以能偵測出極高頻率訊號發生得時間。因 此小波函數實際是一個帶通濾波器，其隨伸縮因數 a 之變動，即對應不同的視窗中心 及視窗半徑。 依小波轉換之定義，可證明連續小波轉換具有以下重要性質[54]： (1) 線性：一個多分量信號之小波轉換等於各分量小波轉換之和。 (2) 自相似性：對應不同尺度參數 a 和不同平移參數 b 之連續小波轉換之間是自相似。 (3) 平移不變性：若𝑓(𝑡)之小波轉換為𝑊𝑓(𝑎, 𝑏)，則𝑓(𝑡 − 𝜏)之小波轉換為𝑊𝑓(𝑎, 𝑏 − 𝜏)。 (4) 伸縮共變性：若𝑓(𝑡)之小波轉換為𝑊𝑓(𝑎, 𝑏)，則𝑓(𝑐𝑡)之小波轉換為𝑊𝑓(𝑐𝑎, 𝑐𝑏)/√𝑐， 其中𝑐 > 0。
17
3.3 建構 ARX 模型
在線性系統中，可將量測自由度間之反應表示成(Huang，1999)[55]： 𝐲(t) = ∑𝐼𝑖=1𝚽𝑖𝒚(𝑡 − 𝑖∆𝑡) + ∑𝐽𝑗=0𝚯𝑗𝒇(𝑡 − 𝑗∆𝑡) (3.13) 其中𝐲(t)為量測動力反應歷時向量（位移、速度、或加速度）；以下假設該向量亦有 n 個分量）；𝒇(𝑡)為系統外力向量歷時；𝚽𝑖與𝚯𝑗待定係數矩陣，∆𝑡為量測反應之取樣時間 間隔；I 及 J 為輸入輸出之延時。以上之離散化運動方程式與多變數 ARX 之時間序列 模型相當類似，式(3.13)為不考慮雜訊或誤差之 ARX 模式。 依 Huang[56]及林[57]之經驗，當量測自由度比系統自由度少很多或有雜訊時，需 用很大之 I 及 J 方能得到較準確之離散化運動方程。但當採用很大之 I 及 J，於利用最 小平方差估算式(3.13)中之Φ𝑖與Θ𝑗時，常會造成數值困難。另外，採用越大之 I 會造成 越多之虛假模態，增加判斷真正系統模態之困難度。若將式(3.13)進行小波轉換，於小 波空間估算Φ𝑖與Θ𝑗；如此便可以適度地降低該數值困難，以較小之 I 及 J 進行估算(Huang and Su，2008)[54]。
將式(3.13)作小波變換，可得： W̅𝒚(𝑎, 𝑏) = ∑𝐼𝑖=1𝚽𝑖W̅𝒚(𝑎, 𝑏 − 𝑖𝑡) + ∑𝐽𝑗=0𝚯𝑗W̅𝒇(𝑎, 𝑏 − 𝑗𝑡) (3.14) 此即小波域中之 ARX 模式。傳統時域下之 ARX 模式僅可藉由調整模型項數( I, J )以 求得最適當之模型；而小波域下之 ARX 模型除了透過模型項數之外，可藉由選擇適 當之小波函數與調整伸縮因數 a 以獲得適當之模型。 在建立式(3.14)時，可利用連續小波轉換之平移不變性；計算𝐲(t)與𝒇(𝑡)之小波轉 換，可同時得𝒚(𝑡 − 𝑖∆𝑡)與𝒇(𝑡 − 𝑗∆𝑡)之小波轉換函數。依小波轉換之特性，於式(3.14) 中取不同之 a 值，代表取某一頻率窗內之反應建立 ARX 模式。為了簡化表示式， W̅𝒚(𝑎, 𝑏 − 𝑖𝑡)與W̅𝒇(𝑎, 𝑏 − 𝑗𝑡)分別記為 𝒚𝑊̅(𝑎, 𝑏̅ − 𝑖)與 𝒇𝑊̅(𝑎, 𝑏̅ − 𝑗)。取特定之 a 與不同 之𝑏̅ 建構式(3.14)，可整理得： [𝑌(0)] = [𝐶] [[𝑌] [𝐹]] (3.15)
18 其中
[𝑌] = [[𝑌(1)]𝑇 [𝑌(2)]𝑇… [𝑌(𝐼)]𝑇]𝑇 (3.15a)
[𝐹] = [[𝐹(1)]𝑇 [𝐹(2)]𝑇… [𝐹(𝐽)]𝑇]𝑇 (3.15b)
[𝑌(𝑖)] = [ 𝒚
𝑊̅(𝑎, max − 𝑖) 𝒚𝑊̅(𝑎, max − 𝑖 + 1) … 𝒚𝑊̅(𝑎, max − 𝑖 + 𝑀)] (3.15c)
[𝐶] = [𝚽1 𝚽2 … 𝚽𝐼 𝚯𝟏 𝚯𝟐… 𝚯𝐽] (3.15d)
max = max( I, J ) (3.15e)
[𝐹(𝑖)] 之定義等同於[𝑌(𝑖)] 欲估算結構系統之動態特性(如：自然振動頻率、模態阻尼比以及振型)，依式(3.13) 知其決定於係數矩陣𝚽𝑖( i = 1,2,…I )。依黃[55]之推導可知，若從Φ 建構矩陣[i 𝐆]： [𝐆] = [ [𝟎] [𝐈] [𝟎] [𝟎] [𝟎] [𝐈] ⋮ ⋮ ⋮ [𝟎] ⋯ [𝟎] [𝟎] ⋯ [𝟎] ⋮ ⋮ ⋮ [𝟎] [𝟎] [𝟎] 𝚽𝐼 𝚽𝐼−1 𝚽𝐼−2 𝚽[𝟎]𝐼−3 ⋯ [𝐈]⋯ 𝚽1] (3.16) 其中，[𝟎]與[𝐈]分別代表維度為N𝑑× N𝑑 之零矩陣與單位矩陣，N𝑑為量測自由度。[𝐆]之 特徵值及特徵向量相關於該結構系統之動態特性。令λ𝑘 = 𝑎𝑘+ 𝑖𝑏𝑘為[𝐆]之第 k 特徵值， 則該結構系統之第 k 模態自然振動頻率 ω𝑘 及模態阻尼比 ξ𝑘 為： ω𝑘 = √𝛼𝑘2+ 𝛽𝑘2 (3.17) ξ𝑘= −𝛼𝑘/𝛽𝑘 (3.18) 其中 𝛼𝑘 =2Δ𝑡1 ln (𝑎𝑘2 + 𝑏𝑘2) ；𝛽𝑘 =Δ 1 tan−1(𝑏𝑎𝑘 𝑘)。
19 令 𝐠𝑘 為[𝐆]之第 k 振態，將 𝐠𝑘 表示成： 𝐠𝑘 = [ {𝜓}1 {𝜓}2 ⋮ {𝜓}𝐼 ] 𝑘 (3.19) 其中 {𝜓}𝑖 含有N𝑑 個分量，由於[𝐆]之特殊結構，可證明： {𝜓}𝑖 = 𝜆𝑘{𝜓}𝑖−1 (3.20) 即不同之 i，{𝜓}𝑖 彼此之間相互平行，故任一{𝜓}𝑖 均對應於結構系統之第 k 模態的振 形。
3.4 以移動最小平方差法架構 TVARX 模型
傳統之基底函數展開法常用之基底函數包含：傅利葉級數(Marmarelis ,1987) [58]、
小波函數(Tsatsanis and Giannakis ,1993)[59]、Legendre polynomials (Niedwiecki ,2000) [60]、Walsh 函數(Zou et. al ,2003)[61]……等。
本節將介紹以移動最小平方差法(MLS)(Liu, 2003)[62]架構 TVARX 模型的方法以 移動最小平方差法(MLS)將 TVARX 之時變係數以基底函數展開之後，先利用權重平 方差估計各個節點(Nodal points)之形狀函數，再以傳統的最小平方差法估算各形狀函 數對應之係數。此方法的優點為能以較少之基底展開，避開用高階多項函數所造成之 數值問題。 由於多項式函數具有平滑性質，對於平滑變化之時變係數的識別結果較佳，因此 本研究乃利用多項式函數，並配合最小移動平方差法，取最高冪次為 2(即[ 1 t t2 ])， 考慮指數權重函數之影響半徑(dm)，求得函數展開之係數。 TVARX 通常被應用於架構時變系統之輸入與輸出之關係，其模型的定義為： 𝐲(t) = ∑𝐼𝑖=1𝚽𝑖(𝑡)𝒚(𝑡 − 𝑖) + ∑𝐽𝑗=0𝚯𝑗(𝑡)𝒇(𝑡 − 𝑗) + 𝒂𝑛(𝑡) (3.21) 其中𝒚(𝑡 − 𝑖)為在t − iΔt 時刻所測量得到之反應向量；𝒇(𝑡 − 𝑗)是在t − jΔt時刻所測量得 之外力輸入向量；𝚽𝑖(𝑡)和𝚯𝑗(𝑡)為待測之時變係數矩陣；𝒂𝑛(𝑡)為由量測噪訊或模型誤
20 差所造成之誤差向量。當利用 TVARX 模型建構非線性系統之輸入及輸出時，即是將 位移或速度反應相關之系統特性(如：阻尼、勁度等)以一等價之時間顯函數來替代。 以移動最小平方差法架構 TVARX 模型，其詳細推導過程可參閱蘇(2008)[54]。在 多自由度時變系統之反應中，TVARX 模型中以多項式基底函數的展開，可以矩陣形 式表示： 𝚽𝑖(𝑡) = ∑𝑁𝑛̂=0𝑖 𝑨𝑖𝑛̂𝑡𝑛̂ = 𝑨̃𝑖𝑷𝑖 (3.22) 𝚯𝑗(𝑡) = ∑𝑁𝑛̂=0̅𝑗 𝑩𝑗𝑛̂𝑡𝑛̂ = 𝑩̃𝑗𝑻𝑗 (3.23) 其中， 𝑷𝑖 = [[𝑰] [𝑰]𝑡 [𝑰]𝑡2 … [𝑰]𝑡𝑁𝑖]𝑇 (3.23a) 𝑻𝑗 = [[𝑰] [𝑰]𝑡 [𝑰]𝑡2 … [𝑰]𝑡𝑁̅𝑗]𝑇 (3.23b) 𝑨̃𝑖 = [ 𝑖0 𝑖1 ⋯ 𝑖𝑁𝑖] (3.23c) 𝑩̃𝑗 = [𝐵𝑗0 𝐵𝑗1 ⋯ 𝐵𝑗𝑁̅𝑗] (3.23d) 𝑨𝑖𝑛̂與𝑩𝑗𝑛̂為待定係數矩陣，[𝑰]為維度與係數矩陣一致之單位矩陣。式(3.22)及(3.23)中
之待定係數𝑨𝑖𝑛̂與𝑩𝑗𝑛̂可以用權重最小平方差法(weighted least square)求解。
令𝚽𝑖(𝑡𝑘)與𝚯𝑗(𝑡𝑘)之真值分別為𝚽̅𝑖𝑘與𝚯̅𝑗𝑘。此外，待定係數可以透過誤差函數最 小化求得，以求取式(3.2a)𝚽中之𝑨̃𝑖為例，將誤差函數定義為： 𝑬̅(𝑡) = ∑ (𝑨̃𝑙̅𝑙=1𝑖 𝑖𝑷𝑖(𝑡𝑙) − 𝚽̅𝑖𝑙)𝑇𝑊(𝑡, 𝑡𝑙)(𝑨̃𝑖𝑷𝑖(𝑡𝑙) − 𝚽̅𝑖𝑙) (3.24) 其中𝑊(𝑡, 𝑡𝑙) = 𝑤(𝑡, 𝑡𝑙)[𝐼]；𝑙̅𝑖為描述𝚽𝑖(𝑡)上所取之節點總數，其值遠小於反應的總取 樣點數。 使函數𝑬̅最小化可得： ∂𝑬̅ ∂𝑨̃𝑖= [0] (3.25) 將誤差函數代入上式，可得：
21 𝑨̃𝑖𝑨̅𝑖(𝑡) = 𝚽̅𝑖𝑸𝑖(𝑡) (3.26) 其中， 𝑨̅𝑖(𝑡) = ∑𝑙̅𝑙=1𝑖 𝑷𝑖(𝑡𝑙)𝑊(𝑡, 𝑡𝑙)𝑷𝑖(𝑡𝑙)𝑇 (3.26a) 𝑸𝑖(𝑡) = [𝑄𝑖1 𝑄𝑖2 … 𝑄𝑖𝑙̅𝑖]𝑇 (3.26b) 𝑸𝑖𝑙 = 𝑷𝑖(𝑡𝑙)𝑊(𝑡𝑛, 𝑡𝑙) (3.26c) 𝚽̅𝑖 = [Φ̅𝑖1 Φ̅𝑖2 ⋯ Φ̅𝑖𝑙̅] (3.26d) 則式中𝑨̃𝑖之解為： 𝑨̃𝑖 = Φ̅𝑖𝑸𝑖(𝑡)𝑨̅𝑖−1(𝑡) (3.27) 上式所得之 ̃𝑖為時間 t 之函數，將之代回(3.26)，可得： 𝚽𝑖(𝑡) = 𝚽̅𝑖𝚽̃𝑝,𝑖(𝑡) (3.28) 𝚽̃𝑝,𝑖(𝑡)類似有限函數中所稱之形狀函數(shape function) 矩陣， 𝚽̃𝑝,𝑖(𝑡) = 𝑸𝑖(𝑡)𝑨̅𝑖−1(𝑡)𝑷𝑖(𝑡)𝑇 (3.29) 同理，𝚯𝑗(𝑡)可用相同方式展開，表示成： 𝚯𝑗(𝑡) = 𝚯̅𝑗𝚯̃𝑝,𝑖(𝑡) (3.30) 其中， 𝚯̅𝑗 = [Θ̅𝑗1 Θ̅𝑗2 ⋯ Θ̅𝑗𝑙̂𝑗] (3.30a) 𝚯̃𝑝,𝑗(𝑡) = 𝑼𝑖(𝑡)𝑩̅𝑗−1(𝑡)𝑻𝑗(𝑡)𝑇 (3.30b)
22 𝑩̅𝑗(𝑡) = ∑𝑙̂𝑙=1𝑗 𝑻𝑗(𝑡𝑙)𝑊(𝑡, 𝑡𝑙)𝑷𝑗(𝑡𝑙)𝑇 (3.30c) 𝑼𝑗(𝑡) = [𝑈𝑖1 𝑈𝑖2 … 𝑈𝑗𝑙̂𝑗] 𝑇 (3.30d) 𝑼𝑗𝑙 = 𝑻𝑗(𝑡𝑙)𝑊(𝑡, 𝑡𝑙) (3.30e) 𝑙̂𝑗為描述𝚯𝑗(𝑡)上所取之節點總數(假設𝛉𝑗(𝑡)中每一分量函數所取之節點總數相同)。在 式(3.28)與式(3.30)中的未知數為𝚽̅𝑖和𝚯̅𝑗，以最小平方差法對𝚽̅𝑖和𝚯̅𝑗求解，其所對應的 誤差函數定義為： 𝐄 = ∑𝑁𝑛̅=1̅ (𝒂𝑛(𝑡𝑛̅))𝑇(𝒂𝑛(𝑡𝑛̅)) (3.31) 其中𝑁̅為用於分析之各自由度反應資料點數，由式(3.21)，(3.28)和(3.30)可得： 𝒂𝑛(𝑡) = 𝒚(𝑡) − [∑𝑖=1𝐼 𝚽̅𝑖(𝑡)𝒚(𝑡 − 𝑖) + ∑𝑗=0𝐽 𝚯̅𝑗(𝑡)𝒇(𝑡 − 𝑗)] = 𝒚(𝑡) − [∑𝑖=1𝐼 𝚽̅𝑖𝚽̃𝑝,𝑖(𝑡)𝒚(𝑡 − 𝑖) + ∑𝐽𝑗=0𝚯̅𝑗𝚯̃𝑝,𝑗(𝑡)𝒇(𝑡 − 𝑗)] = 𝒚(𝑡) − [∑𝑖=1𝐼 𝚽̅𝑖𝚪𝑡,𝑖+ ∑𝐽𝑗=0𝚯̅𝑗(𝑡)𝛀𝑡,𝑗] (3.32) 其中， 𝚪𝑡,𝑖 = [𝚽̃𝑝,𝑖,1𝒚(𝑡 − 𝑖)𝑇 𝚽̃𝑝,𝑖,2𝒚(𝑡 − 𝑖)𝑇 ⋯ 𝚽̃𝑝,𝑖,𝑙̅𝑖𝒚(𝑡 − 𝑖) 𝑇]𝑇 (3.32a) 𝛀𝑡,𝑗= [𝚯̃𝑝,𝑗,1𝒚(𝑡 − 𝑗)𝑇 𝚯̃𝑝,𝑗,2𝒚(𝑡 − 𝑗)𝑇 ⋯ 𝚯̃𝑝,𝑗,𝑙̅𝑗𝒚(𝑡 − 𝑗)𝑇]𝑇 (3.32b) 對未知數取導數可得： ∂𝐄 ∂𝚽̅𝑖 = ∑ 2𝚽̃𝑖𝒚(𝑡𝑛̅− 𝑖)𝑇 𝑁̅ 𝑛̅=1 [𝒚(𝑡𝑛̅) − (∑ 𝚽̅𝑖𝚪𝑡,𝑖+ ∑ 𝚯̅𝑗(𝑡)𝛀𝑡,𝑗 𝐽 𝑗=0 𝐼 𝑖=1 )] = [0]
23 (3.33a) ∂𝐄 ∂𝚯̅𝑖 = ∑ 2𝚯̃𝑖𝒚(𝑡𝑛̅− 𝑗)𝑇 𝑁̅ 𝑛̅=1 [𝒚(𝑡𝑛̅) − (∑ 𝚽̅𝑖𝚪𝑡,𝑖 + ∑ 𝚯̅𝑗(𝑡)𝛀𝑡,𝑗 𝐽 𝑗=0 𝐼 𝑖=1 )] = [0] (3.33b) 將以上聯立方程組以矩陣形式表示之，可得： 𝑽𝑚𝑇𝑭 = 𝑽𝑚𝑇𝑽𝑚𝑪𝑚 (3.34) 其中， 𝑽𝑚 = [ 𝚪𝑡1,1 𝚪𝑡1,2 … 𝚪𝑡2,1 𝚪𝑡2,2 … ⋮ ⋮ ⋱ 𝚪𝑡1,𝐼 𝛀𝑡1,0 𝛀𝑡1,1 𝚪𝑡2,𝐼 𝛀𝑡2,0 𝛀𝑡2,1 ⋮ ⋮ ⋮ … 𝛀𝑡1,𝐽 … 𝛀𝑡2,𝐽 ⋱ ⋮ 𝚪𝑡𝑛̅,1 𝚪𝑡𝑛̅,2 … 𝚪𝑡𝑛̅,𝐼 𝛀𝑡̅𝑛,0 𝛀𝑡̅𝑛,1 … 𝛀𝑡̅𝑛,𝐽 ] (3.34a) 𝐅 = [𝒚(𝑡1) 𝒚(𝑡2) ⋯ 𝒚(𝑡𝑛̅)]𝑇 (3.34b) 𝑪𝑚 = [𝚽̅1𝑇 𝚽̅ 2𝑇⋯ 𝚽̅𝐼𝑇 𝚽̅0𝑇 𝚽̅1𝑇 ⋯ 𝚽̅𝐽𝑇] 𝑇 (3.34c) 對式(3.34)求解即可得到每個節點之真值，代回式(3.28)即可得到此系統於某個時刻 t=𝑡𝑛之時變係數𝚽𝑖(𝑡𝑛)。
3.5 瞬時模態參數估算
架構 TVARX 之各個時變係數時，若取相同之多項式階數以及相同之節點數(即令 𝑁1 = 𝑁2 = ⋯ = 𝑁𝐼 = 𝑁̅0 = 𝑁̅1 = ⋯ = 𝑁̅𝐽 = 𝑁𝑝；𝑙1 = 𝑙2 = ⋯ = 𝑙𝐼 = 𝑙̂0 = 𝑙̂1 = ⋯ = 𝑙̂𝐽 = 𝑙𝑛)，則每個時變係數可以使用相同之形狀函數，使得運算過程更有效率。 TVARX 通式如式(3.21)所示，在瞬時 t 下，Φ𝑖(𝑡)與Θ𝑗(𝑡)均為常數矩陣，故 TVARX 模式與非時變之 ARX 模式對等。 依非時變 ARX 模式估算動態特性之方法，令24 [𝐆] = [ 0 𝑰 0 0 0 𝑰 ⋮ ⋮ ⋮ 0 ⋯ 0 0 ⋯ 0 ⋮ ⋱ ⋮ 0 0 0 𝚽𝑖(𝑡) 𝚽𝑖−1(𝑡) 𝚽𝑖−2(𝑡) 𝚽𝑖−30 (𝑡) ⋯ ⋯ 𝚽𝑰 𝐼(𝑡)] (3.35) 時變系統中，[G]隨時間而改變。[G]之特徵值及特徵向量式與結構系統之瞬時動態特 性有關，同非時變系統之推導可知： [𝐆]𝛙̃𝑘= 𝜆𝑘𝛙̃𝑘 (3.36) 其中𝜆𝑘及𝛙̃𝑘為[G]之第 k 個特徵值及特徵向量。 在時變系統中，𝜆𝑘及𝛙̃𝑘亦是時間函數，令： 𝛙̃𝑘 = { 𝛙̃𝑘 (1) 𝛙̃𝑘(2) ⋮ 𝛙̃𝑘(𝐼)} (3.37) 其中，𝛙̃𝑘(𝑖)為一(𝑙 ×1)之向量，l 為量測自由度，且由[G]可知： 𝛙̃𝑘(𝑖) = 𝜆𝑘𝛙̃𝑘(𝑖−1) (i = 1 ,2 , …,I ) (3.37a) ∑𝐼−1𝑗=0𝝋̃𝑛−𝑗𝛙̃(𝑗+1)𝑘 = 𝜆𝑘𝛙̃𝑘(𝐼) (3.37b) 由(3.37a)及(3.37b)可得： 𝛙̃𝑘= { 𝛙̃𝑘(1) 𝜆𝑘𝛙̃𝑘(1) 𝜆𝑘2𝛙̃𝑘(1) ⋮ 𝜆𝑘𝐼−1𝛙̃𝑘(𝐼)} (3.38) 及𝛙̃𝑘(1)代表量測自由度之第 k 瞬時模態。 式(3.35)中之特徵值常為複數或成對之共軛根。令𝜆𝑘= 𝑎𝑘+ 𝑖𝑏𝑘，則結構系統之瞬
25 時擬自然振動頻率及阻尼比為： 𝛽̃𝑘 = √𝛼𝑘2+ 𝛽𝑘2 ， 𝜉𝑘 = 𝛼𝑘 𝛽̃𝑘 ⁄ (3.39) 其中，𝛽̅𝑘 = Δ𝑡1 𝑡𝑎𝑛−1 𝑏𝑘 𝑎𝑘 ， 𝛼𝑘 = 1 2Δ𝑡ln (𝑎𝑘 2 + 𝑏 𝑘2) (3.40) Δt為時間增量(及取樣頻率之倒數)，通常𝛽̃𝑘亦稱為瞬時擬自然振動頻率( instantaneous
26
第四章 數值模擬之系統識別結果分析
4.1 前言
本研究的數值模擬結構系統受地震反應分為兩大部分討論，分別是非時變系統( Time-invariant system )與時變系統( Time-variant system)考慮濾波對於系統識別的影響。 由於時變系統的行為較非時變系統複雜得多，因此如果濾波方法及系統識別適用於時 變系統的話，那麼該分析方法同樣也應適用於非時變系統；可知若是分析方法對於非 時變系統是可行的，則我們將可以試著將其方法套用至時變系統，看看是否同樣可行； 但假若該分析方法不適用於非時變系統，那麼將同樣也不適用於時變系統，分析流程 如圖 4.1 所示。 為瞭解本研究所採用之濾波方法對於非時變系統與時變系統的識別效果，將透過 一系列之數值模擬反應對此流程進行驗證：非時變系統使用七層樓之剪力結構物之數 值模擬地震反應，建構 ARX 模型以進行系統識別；時變系統使用五層樓剪力構架之 數值模擬地震反應建構 TVARX 模型；接著再探討信號經過濾波之後的識別結果與未 濾波前的識別結果之差異。4.2 節分別介紹非時變系統的無噪訊反應(NSR=0%，NSR 為噪訊比，即 Noise-Signal Ratio)和含噪訊反應(NSR=10%、20%、30%、40%)的識別 結果；4.3 節介紹非時變系統經過第二章所述之三種濾波方法濾波過後的系統識別結果； 4.4 節分別介紹時變系統的無噪訊反應和含噪訊反應(NSR=1%、3%、5%)以及其經過 濾波器處理後的識別結果。
4.2 非時變系統之數值模擬的系統識別分析結果
圖 4.2a 與圖 4.2b 為本研究中七層樓之剪力構架示意圖，其數值模擬地震反應的取 樣頻率均為 250Hz，且每層樓的質量 mi為 0.1 噸，勁度 ki為 80kN/m，其中 i=1~7。圖 4.3a~圖 4.3g 所示為非時變系統之七層樓剪力構架之第一層樓到第七層樓中，其訊號不 含噪訊與含不同噪訊比時之反應。4.2.1 小波函數之影響
本節參考蘇(2008)[54]，介紹以 Meyer 小波函數為基底來架構 ARX 模型時，a 值 的選取方法。Meyer 小波函數具有頻率轉換函數之數學表示式，且具有很類似帶通濾
27 波器之特性。因此，以下探討利用 Meyer 函數處理含有噪訊的反應。 Meyer 小波之小波函數M定義於頻率域： Ψ̂𝑀(𝜔) = { (2𝜋)−1/2𝑒𝑖𝜔/2sin [𝜋2𝜈 (2𝜋3 |𝜔| − 1)] ，2𝜋3 ≤ |𝜔| ≤4𝜋3 (2𝜋)−1/2𝑒𝑖𝜔/2cos [𝜋 2𝜈 ( 3 4𝜋|𝜔| − 1)] ， 4𝜋 3 ≤ |𝜔| ≤ 8𝜋 3 0 ， |𝜔| ≤ [2𝜋3 ，8𝜋3] (4.1) 其中，ν(s)為構造 Meyer 小波之輔助函數： ν(s) = 𝑠4(35 − 84𝑠 + 70𝑠2− 20𝑠3) 𝑠 ∈ [0,1] (4.2) Ψ̂𝑀(𝜔)為有限支撐長度之正交小波。因此，配合此小波函數之”主要頻率保留區間”， 即可決定出適當之 a 值。 由式(4.1)顯示，Ψ̂𝑀(𝜔)在2 38 3之區間才有值。若定義頻率分量為最大 者 之 90% 以 上 之 頻 率 區 間 為 ” 主 要 頻 率 保 留 區 間 ” ， 由 式 (4.1) 可 得 此 區 間 為 [0.5348,0.9311]。因此，對尺度因數 a 之 Merey 小波函數的”主要頻率保留區間”範圍為 [0.5348/a,0.9311/a]。此小波函數之頻率域中，除了定義之區間外，其餘部份均為 0。 由此可知，除非選擇 a 值所對應之”主要頻率保留區間”涵蓋欲識別之振動頻率，否則 將難以獲取準確之結果。
