結構健康診斷及控制研究---大型結構實驗驗証---子計畫：應用類神經網路技術於房屋結構之非線性系統識別、損壞診斷以及健康監測(II)

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

■ 成 果 報 告

□期中進度報告

結構健康診斷及控制研究:大型結構實驗驗証-子計畫:應用類神經網路技術於

房屋結構之非線性系統識別、損壞診斷以及健康監測(II)

計畫類別：□ 個別型計畫 ■ 整合型計畫

計畫編號：NSC

94-2625-Z-009-005-執行期間： 94 年 8 月 1 日至 95 年 7 月 31 日

計畫主持人：洪士林

共同主持人：

計畫參與人員： 林宏宇、張巍贏

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

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

□赴國外出差或研習心得報告一份

□赴大陸地區出差或研習心得報告一份

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

□國際合作研究計畫國外研究報告書一份

處理方式：除產學合作研究計畫、提升產業技術及人才培育研究計畫、

列管計畫及下列情形者外，得立即公開查詢

□涉及專利或其他智慧財產權，□一年□二年後可公開查詢

執行單位：國立交通大學 土木工程學系(所)

中 華 民 國 95 年 12 月 1 日

摘 要

結構物的系統參數如：自然頻率（Natural Frequency）、振態阻尼比（Damping Ratio） 為結構物的重要特性之一，若能得到結構物的系統參數，將能掌握結構物在風力或地震力 影響之下的動力行為表現，亦可運用在結構物的健康檢測。本研究主要是採用倒傳遞類神 經網路（Back-Propagation Neural Network, BPN）及希伯特-黃轉換兩種方法來識別結構 物參數。倒傳遞類神經網路進行結構系統模態識別的架構主要為類神經網路及時間域分析 技術所組成。使用地表加速度的資料以及各樓層的加速度反應作為類神經網路的訓練案 例，設定與調整網路之適當的層數、自由度等參數後，配合著類神經網路的學習能力與疊 代計算可得到網路權值。將訓練完畢的網路權植，代入動力方程式中即可估算出系統之動 力特性。由黃鍔等人所提出的希伯特-黃轉換(Hilbert Huang transform,HHT)，包含了經 驗模態分解法（Empirical Mode Decomposition,EMD）與希伯特轉換（Hilbert

Transform,HT）。經驗模態分解法將各感測器接收到的真實地震訊號分解成數個重要的內建

模態函數（Intrinsic Mode Functions，IMF），配合使用帶通濾波器可得到各個模態的資

訊，再由隨機遞減法（Random Decrement Technique,RDT）將得到的各個模態的反應訊號 轉換成為自由震動衰減的訊號，最後以希伯特轉換得到瞬時的相位角與振幅，經由動力方 程式的推導可求出結構物的阻尼比與自然頻率。本研究先以建立好的數值模型來測試此兩 種方法，驗證是否皆能得到預期的答案，隨後以Benchmark模型在振動台上試驗的反應資料 來進行測試，比較兩者所求出真實結構物的系統參數是否相近。在自然頻率識別的部分， 兩者皆有不錯的是識別效果，而在阻尼比的識別上類神經網路依然可以得到準確的答案， 希伯特黃轉換則會因為參數的改變得到差異較大的結果。然而，類神經網路對於雜訊過高 的實驗案例會出現無法識別的現象，而希伯特黃轉換可藉由調整濾波器的參數得到識別結 果。 關鍵詞: 類神經網路；希伯特黃轉換；經驗模態分解法；系統識別；基準模型

Abstract

The modal parameters of a structure, such as natural frequency and damping ratio, are the important characteristics of the structure as they are the key factors for the technology of structural health moniroing. Herein, two different approaches, artificial neural networks (ANNs) and Hilbert-Huang transform (HHT), are employed to identify these parameters of the structure under different external excitations.

The process of system identification using back-propagation neural network (BPN) consists of two submodels, artificial neural network model and time domain analysis technique. First, the dynamic responses of each floor and external excitation are used to train the ANN model. Following, the natural frequency and damping ratio of the structure can be computed via the time domain analysis with the trained weights of the network.

Hilbert-Huang transform consists of Empirical Mode Decomposition (EMD) and Hilbert Transform (HT). EMD is the method of signal processing, it can decompose the inputted signal to several Intrinsic Mode Functions (IMFs). Passing the first IMF to the bandpass filter, each mode responses can be obtained. The Random Decrement Technique is used to get the free vibration signal. Afterwards, the instantaneous phase angle and amplitude can be got by passing free vibration signal through HT. Following, the natural frequency and damping ratio of the structure can be computed via the derivation of dynamic equations.

In this study, a three-floor benchmark model with different configuations is employed as a testing structure. The numerical analysis and corresponding shake table test results with different configurations of the benchmark model are used to verify these two methods and the verification results are also compared. The results indicated that both approaches have good performance for identifying nature frequencies of variant benechmark models. For identification of damping ratios, ANN model also performs well, but HHT does not.

Keywords: Artificial neural networks (ANNs); Hilbert Huang transform (HHT)；Empirical Mode Decomposition (EMD; System identification (SI); Benchmark model

壹、前 言

二十一世紀的今天，隨著商業與科技的蓬勃發展，人類文明欲趨都市化，而產生顯著 的人口集中現象，人們為了追求更方便與有效率的空間使用，以及工程師對高樓建築極限 的挑戰，致使都會城市中高樓大廈林立，且伴隨著各大都會區之間的人車頻繁往來，城鄉 機場港口之間的人力物流運輸，長距離高流量的交通需求大增，也使得捷運、大型橋樑、 高速鐵路...等，大量的興建。而台灣正好處於環太平洋的地震帶西緣，因菲律賓板塊與歐 亞大陸板塊長期處於碰撞狀態下，平均每年發生百餘次的有感地震。若在都市人口密集的 地區，發生高規模的大地震，造成的傷亡與財物損失是令人難以想像的，幾年前的921大地 震就是一個令人難以抹滅的慘痛經驗。為了避免地震對人們傷害，對結構物做非破壞檢測 以進行補強修復，選擇使用合適的隔減震設備也將是勢在必行的；而在這些措施的執行使 用之前，我們必須先對結構物做系統識別，來確認結構物的系統參數，如此才能對地震做 出最有效的防範。 在本計畫前一階段的報告中，有多筆 Benchmark 模型在振動台上所量測的地震反應， 以這些數據做為系統識別資料，應用兩種發展出來的方法，來比較是否能正確得到相同的 系統參數。

1.1、研究背景

在過去的二十年間，類神經網路（Artificial Neural Network, ANN）已經逐漸發展 成熟並且成功的應用在影像、文字、語音的辨識系統、礦床探測、製程監控、股票預測、 財務分析、氣象預測、疾病診斷、土木工程…等領域之中。類神經網路主要具有下列之優 點： 1. 具有過濾能力與高容錯性： 在類神經網路中各個輸入與輸出的關係，並不是直接由網路中單一個節點所負責的，而 是每一個節點都會映射給輸入-輸出模式一組權值。當網路中所有結點與對應的權值組 合在一起時，才能夠映射出完整的輸入－輸出模式。因此，當某一節點所輸入要處理的 訊號具有雜訊或處理的數據不完整時，此一輸入對網路所造成的影響，將不會導致嚴重 的錯誤。 2. 具有適應性學習能力： 每種類型的類神經網路都有其特定的學習演算法，經由演算法可調整節點與節點之間的 連接權重值，藉由不斷地調整權重值，直到產生最佳的輸入－輸出模式，這種能力稱之 為適應性學習能力。 3. 多輸入輸出之系統： 輸入層與輸出層可以具有任意數目的節點，可方便我們運用在多自由度的系統上。 因 此使其無論在訊號處理或系統識別控制及複雜的映射（mapping）等問題上，成為了一 種強有力的工具。

(NASA) 大華府地區學者黃鍔博士近幾年於英國皇家學會(The Royal Society)期刊上所發 表的發明[1]，該研究榮獲 NASA 所頒發的年度發明獎，也因該研究的成就，使其於 2000 年獲選為美國國家工程院院士，2004 年獲選為台灣中央研究院院士。希伯特黃轉換的優勢 在於可處理非線性和非穩態資料，過程主要包含了兩個部分，有前處理的經驗模態分解法 (Empirical Mode Decomposition, EMD)以及之後的希伯特轉換（Hilbert Transform, HT）。 此方法目前同樣已經廣泛的運用於氣候、水利、地震工程、地球物理探測、潛艇設計、結 構損害偵測、衛星資料分析、血壓變化和心律不整等各項研究之中。 因此本研究的主要目的在於建立以此兩種方法為基礎的結構物系統識別程序，並驗證 比較這兩種方法的精確性，且於實地案例中能否得到同樣的結果。

1.2、文獻回顧

如果把結構物視作一個系統，則可利用結構矩陣與結構動力等概念建立出其等效的數 學模型，當此系統接收到一個外部的輸入訊號時，如地震力、風力…等，將可從結構物的 實際量測或數值模型的推導對應得到此一結構物反應的輸出訊號，如：樓層的位移、速度、 加速度。而所謂的系統識別就是利用輸入訊號與輸出訊號之間的關係反推出其等效的數學 模型進而得到此結構物系統的各項參數。所求得的參數可用來檢驗結構的建造是否如同我 們設計時所預期，亦可以用來評估建築物經歷過地震、強風或材料的老化後的強度及其安 全性，可做為補強的依據，來確保結構物的安全。 一般的結構系統識別依據描述系統的方式，可分為參數模式[2]與非參數模式識別。若 從識別時所用的資料，可分為時間域與頻率域識別。又依量測數據之處理程序，則系統識 別技巧可再概分成兩類，線上模式（on-line model）及離線模式（off-line model）。所 謂參數識別是指識別代表該結構物系統的數學模式中，具有物理意義之參數，像是一般在 結構動力分析裡的自然頻率ω及阻尼比ξ等。至於非參數的識別乃對系統不作任何物理或 數學模型之假設，僅對已知的輸入和輸出資料之間以一組相關的數學函數來連結，而此函 數的係數，即為識別的對象，是沒有明確的物理意義。時間域的系統識別是直接將所量測 的時間域反應來做系統識別的資料，而頻率域系統識別法則是先將時間域反應的資料轉換 到頻率域再進行系統識別，頻率域識別若在時變或非線性系統，可以移動視窗方式找尋等 值線性系統參數，優點為計算快速，缺點則是精確度受限於所收集資料長短和取樣頻率大 小。時間域識別，需選擇適當數學模式和合理的誤差準則(Error Criterion)，再運用數值 方法識別出相關參數。線上模式為數據進來一筆；即進行系統識別一次，其分析技巧架構 於遞迴公式（Recursive Formula）[3]。而離線模式，則是將所有量測到的數據記錄完成 後，再進行系統識別。線上分析模式，雖有減少數據儲存時間、計算快速及可分析時變系 統之優點，但於分析線性非時變系統，則其精度一般不如離線分析模式。本研究所應用的 兩種方法，都是屬於離線的時間域參數識別。 許多的系統識別方法已成功地應用於土木工程，由於估算方法之多樣性，有不少專家 學者即針對各種方法進行比較，如Yun 和Shinozuka [4]則針對傳統最小平方差法、工具變 數法、最大可能機率法(Maximum Likelihood Algorithm)和有限訊息最大可能機率法

(Limited Information Maximum Likelihood Method)應用於多重輸入╱輸出系統建立 ARX(Auto Regressive eXogenous) 模式之比較。Ghanem 和Shinozuka [5]、[6]亦對推廣 卡氏過濾法(Entended Kalmen Filter Algorithm)、最大可能機率法、遞迴最小平方差法 (Recursive Least Square Method)、和遞迴工具變數法應用於多重輸入╱輸出系統建立ARX 模式比較。Saridis [7]針對線上分析法之互相關法（Cross Correlation Algorithm）、 一階隨機近似法（First Order Stochastic Approximation Algorithm）、二階隨機近似 法(Second Order Stochastic Approximation Algorithm)、最大可能機率法、最大事後機 率法（Maximum Posteriori Probability Algorithm）、推廣卡氏過濾法應用於線性單一 輸入╱輸出系統(SISO)之比較。

類神經網路之相關研究已有近五十年的歷史，最近二十年逐漸發展成熟並應用在許多 的領域之中。Rosenblatt（1957）提出第一種神經網路模式-感知機（Perceptron）模式， 它由二元值神經元組成，以此模仿生物的大腦及視覺系統，主要用於理論研究與樣本識別。 Widrow和Hoff（1960）提出自適應線性元件（Adaptive Linear Element）模式，它是一種 連續值的線性網路。Rumelhart等人（1985）提出著名的倒傳遞類神經網路模式

（Back-Propagation Network, BPN），此網路模式之基本原理是利用最陡坡降法（Gradient Steepest Descent Method）的觀念，將誤差函數予以最小化。Kohonen（1980）提出組織 映射圖（Self-Organizing Map）模式，並且在1988從該模式衍生出學習向量量化網路 （Learning Vector Quantization Network）。Hinton 和Miller（1988）對霍普菲爾-坦 克神經網路在解決最佳化問題時，其收斂最小值和參數設定等問題加以改善。Barnard （1992）探討訓練類神經，網路值目標函數最佳化之各種方法，並提出一個以隨機的觀念 所建立的序列演算法。Hagan和Menhaj（1994）根據Kollias和Anastassiouh所建議的概念， 提出改良式倒傳遞演算法，此法是將應用與非線性最小平方法的Levenberg-Marquardt演算 法， 配合傳統倒傳遞演算法來訓練前向式類神經網路。Flood和Kartam（1994）提出將類 神經網路應用於土木工程上的明瞭性、使用性以及實用性方面的論述，並利用倒傳遞前向 式網路訓練，來解決結構分析問題。Narendra和Parthasarathg(1990)曾驗證類神經網路能 有效地使用至非線性動力系統識別中。 Wu 等人[8]以一系列三層樓鋼構架之數值模擬資 料，利用倒傳遞神經網路（Back-Propagation Neural Network, BPN）來描述該結構的破 壞狀態。該研究以加速度反應富氏譜以及桿件勁度，分別作為其BPN 之輸入及輸出變數。 而Elkordy 等人[9]則以模態作為其BPN 之輸入變數，以偵測模擬的結構破壞。Szewczyk 與 Hajela [10]則應用反傳遞神經網路（Counter-Propagation Neural Network, CPN），以 剛架之靜定位移來估算桿件勁度的折減情形。Pandey 與Barai [11]應用多層感知器 （Multilayer Perceptron），以數值模擬資料偵測桁架橋之破壞。Zhao 等人[12]以靜定 位移、自然頻率、以及模態，應用CPN 來分別偵測梁和剛架的破壞位置。Masri 等人[13] 根據非線性系統識別，建立了一套破壞偵測的方法。其方法裡採用了實驗中所量測到的位 移、速度、加速度反應，以及輸入外力等資料，作為網路訓練之用。傳統的類神經網路多 為一非參數系統識別之方法，而本文則是利用ARX來模擬結構系統的運動方程式，建立出參 數系統識別之方法，而推出結構的系統參數。

線性且穩態訊號的技術，及之後發展的像是小波和Wanger-Ville distribution [14][15] 用來處理非穩態但線性的資料，還有像是Tong [16]、Kantz和Schreiber [17]、Diks [18] 提出數種可以分析非線性但穩態的系統，希伯特黃轉換則是一種可以處理非線性非穩態時 間序列資料的分析方法。其主要包含了經驗模態分解法與希伯特轉換兩個步驟。EMD的過程 是利用資料變化的內部時間為尺度，將資料分解成多個內建模態函數(Intrinsic Mode Function, IMF)的疊加，這樣可以看作是將原來的訊號以IMF 當做基底展開，它也可以是 非線性或非平穩性的，完全基於原來訊號的性質。這樣的過程可以使原始訊號的複合波拆 解成單一波，如此一來對每個IMF都將具有一個良好的希伯特轉換對，而可得到訊號的頻 率、振幅、相位角。Ruzzene [19]等證明當模態的自然頻率非常接近時，希伯特黃轉換比 小波更有能力處理，Yang和Lei [20]提出了一種以HHT為基礎的方法來分析多自由度的線性 結構，Vincent、Hu和Hou [21]比較經驗模態分解法與小波分析在結構物的損害偵測，Yang 和Lei [22]～[25]延伸先前的研究，利用自由振動的加速度資料，不止能識別線性結構物 的自然頻率、阻尼係數更能推演到質量、勁度、阻尼矩陣。Yang和Lei [26]提出了整合使 用經驗模態分解法、隨機遞減法、希伯特轉換處理高樓風力的振動資料，來做系統識別， Yang、Lin、Pan和Huang [27][28]等人將經驗模態分解法搭配帶通濾波器來求得各模態的 系統參數，Yang、Lin與Pan [29]應用希伯特黃頻譜來檢測的損害程度，XU、CHEN和ZHANG [30] 使用希伯特黃轉換搭配隨機遞減法分析深圳的地王商業中心於強颱約克的侵襲下所記錄的 資料，而得到結構物的系統參數，Yang、Lin和Huang [31]運用希伯特黃轉換來識別結構物 損傷後的自然頻率和阻尼比，Banfu Yan 和 Ayaho Miyamoto [32]比較小波與希伯特黃轉 換在系統識別上優缺點。Huang & Shen [33]又提出了完整的希伯特黃轉換原理解說與應用 及其還之後要深入研究改進的部分。而針對實際的地震資料目前尚未有人提出具體的測試 結果，本研將以Benchmark振動台的資料為測試的案例，來驗證此方法的可行性。

貳、類神經網路於系統識別之理論與架構

類神經網路（ANN）近年來被快速地廣泛發展應用於各領域，主要其乃是一個具有計算 能力的多層網路系統，原理在於使用大量簡單而具有平行處理能力的人工神經元來模擬人 腦的學習行為，所以常被拿來應用於人工智慧上，而本文將利用其學習能力，在多次的疊 代後，識別出結構物的系統參數。本章首先是倒傳遞類神經網路的介紹，接著是 L-BFGS 演 算法與線搜尋演算法之理論分析，最後是運用於系統識別之架構。

2.1、倒傳遞類神經網路

類神經網路是由生物神經網路得到靈感所發展出的一種系統；它由一些相互連結在一 起的簡單處理單元（結點）所組成，連結的權值（weights）代表儲存在系統的資訊並用來 表示連結的強度，這些權值掌握了使類神經網路產生功能的關鍵。在各種不同的類神經網 路模式中，使用誤差倒傳遞演算法之向前饋入、多層、監督式的神經網路，即所謂的倒傳 遞網路（Back-Propagation Network, BPN）[34]，由於它的簡單性，是目前應用最普遍的 類神經網路學習模式。一個類神經網路在可以應用之前，必須先從一筆存在有一對輸入值 及輸出值的案例或資料庫來訓練。 如圖(2.1)所示，BPN 的網路架構包含了一層的輸入層、一或多層的隱藏層以及一層的 輸出層。而每一層之節點皆與其鄰層的節點相連接。通常隱藏層之結點數目越多收斂越慢， 但可達到更小的系統誤差值，當超過一定數目後，再增加則對降低系統誤差幾乎沒有幫助， 只是徒然增加執行之時間。另外值得一提的是，Hecht-Nielsen [35]在其研究中證明，一 層的隱藏層已足夠解決大部分實際應用上的問題。因此，於本報告中之各個神經網路將只 使用一層的隱藏層。在一類神經網路能夠使用之前，它必須先經過訓練的過程。利用 BP 學 習演算法的訓練過程，一般包含了三個階段。第一階段稱之為資料向前饋入（data feedforward）。輸出層中第i個節點的計算輸出值yi定義如下，見圖（2.1）： 1 1 ( h ( ( i ) )) N N i ij jk k j wi j k y g w g v x

θ

_ν

θ

= = =

∑

∑

+ + , i=1, 2, _L N_o, (2.1) 其中 wij為隱藏層及輸出層節點之間的連接權植；vjk為輸入層及隱藏層節點之間的連接權 值；

θ

_wi與

θ

_vj為轉換函數 g 之門檻值；xk為輸入層第 k 個節點的輸入值。而 Ni、Nh、及 No 則分別為輸入層、隱藏層、及輸出層的節點數目。而轉換函數之採用可取線性或非線性。 第二階段稱之為誤差向後推導（error back-propagation）。在訓練的過程中，以一系 統誤差函數來監測網路的學習表現。而此函數通常定義如下： 1 1 ( ) ( )( ) 2 P T p p p p p E P = =

∑

− − W Y% Y Y% Y , (2.2) 其中P為學習的案例數。 o 1 2 ( y y yi y )N = Y% % % L L% % ； ( _{1 2} y ) o i N y y y = Y _{L L} ，

y

%

_i為輸出節點i

之期望值，而 11 12 1 2 11 12

(

h i h jk N N N ij

v v

v

v

θ θ

_{ν ν}

θ

_ν

w w

w

=

W

_{L L}

_L

_{L L}

_{1 2} ) o h o N N w w wN w θ θ _Lθ 。 訓練的最後階段為權值的修正。標準 BP 演算法係基於最陡梯度法（gradient descent method）並使用固定的搜尋步幅（step length）或學習速率（learning ratio） 來訓練 網路。其權值的修正如下： (k+1) = ( )k + Δ ( )k W W W (2.3) ( ) ( ) k k E η ∂ Δ = − ∂ W W (2.4) 其中η為學習率，一般介於 0~1 之間。上標(k)表示迭代第k次，亦即網路經過k次的學習。 BP 演算法中最小化搜尋方向是由負的誤差函數梯度決定，這種搜尋方向上的搜尋步幅由固 定的學習速率決定，因此常常導致學習之系統誤差不穩定以及學習速度緩慢之困擾。

2.2、L-BFGS 類神經網路

由於 BP 神經網路通常需要大量的學習時間，而且其網路收斂速度非常倚靠學習率的選 擇 。 因 此 ， 本 報 告 中 ， 將 採 用 另 外 一 種 基 於 最 小 記 憶 體 之 BFGS （Broyden-Fletcher-Goldfarb-Shanno）[36]假牛頓二階（Quasi Newton second-order） 法，且結合非精確的線搜尋演算法（inexact line search algorithm），所發展的可調式

L-BFGS 學習演算法[37]，來使得學習的過程更具有效率。在傳統的 BFGS 法中，誤差函數E 的 Hessian 反矩陣 Hk+1可經由下式而逼近。 1

(

)

(

)

T T T k k k k k k k k k k k T T k k k k k k +

= −

−

+

≡

+

H

I ρ s y H I ρ y s

ρ s s

V H V

ρ s s

(2.5) 其中

1/

T k

=

k k

ρ

y s

(2.6) T k

= −

k k k

V

I ρ y s

(2.7) 1 k

=

k+

−

k

s

W

W

(2.8) 1 k

=

k+

−

k

y

g

g

(2.9) k

E

∂

=

∂

g

W

(2.10) 在本法中，僅需儲存 sk及 yk向量，而不像在 BFGS 方法中需建立矩陣 Hk。首先先將這些向 量定義好，然後從最近幾次迭代的資料中，動態地（dynamically）更新近似 Hessian 矩陣

之反矩陣。其優點為簡化計算與減少儲存空間。因此，最後一個步驟之權值調整修正如下： (k 1) ( )k k k

α

+

₌

₊

W

W

d

(2.11) 其中α_k為搜尋步幅；而 dk為搜尋方向，其定義如下 1 k

= −

k k

+

β

k k−

d

H g

d

(2.12) 其中 ( 1) ( 1) ( 1) ( 1) ( 1) T k k k k T k k

β

− − − − −

=

y

H

g

y

d

(2.13) L-BFGS 學習演算法在學習的過程中非精確線搜尋演算法來調整搜尋步幅，而不使用固定之 學習速率。非精確線搜尋演算法是基於下列三個連續的步驟：建立區間（Bracketing），分 割（Sectioning），及內插（Interpolation）。因此，搜尋步幅α_k在每一次迭代過程中需 要滿足下列條件[37]： ( ) ( ) ( ) ( ) ( ) ( ( ) ) (0,1) 0 k k k T k k k k k E E E and

α

βα

β

α

+ ≤ + ∇ ∈ > W d W W d (2.14) ( ) ( ) ( k )T ( ( k )T ( ,1) 0 k k k k k E

α

θ

E

θ

β

and

α

∇ W + d d ≥ ∇ W d ∈ > (2.15) ( ) ( 1) ( k )T 0 k k k E

α

₊ ∇ W + d d < (2.16) 所以，可調式的 L-BFGS 學習演算法可解決標準 BP 演算法裡以試誤法選擇學習速率的問題。

2.3、線性結構系統之類神經網路識別架構

在一般之地震反應監測系統中，通常所量測的多為結構加速度反應。因此，以這些量 測的加速度資料用來訓練類神經網路，待該網路測試收斂至穩定時，將其學習完畢之權值， 拿到系統識別的程序之中，便可輕易的得到系統的動態特性。 本研究所採用的類神經網路架構如圖(2.1)所示。其中 ) (t i f_l − _, i = 0 ,1,2,_Ln _(2.17) 表示在(t-i)時刻下，第l個輸入自由度之地表加速度。而 ) (t j x&&_k − _, j =0 ,1,2,_Lm _(2.18) 表示在(t-j)時刻下，第k個自由度的加速度觀測量。因此，由網路架構圖可知，本報告中 採用了三層的倒傳遞網路架構，輸入層取地震前(t-i)個時刻的外力輸入與前(t-j)個時刻 的觀測自由度加速度反應值，作為網路輸入變數，輸出層則取 t 時刻各觀測自由度的加速度 反應值，作為網路的輸出變數，而隱藏層為一層。學習演算法則採用了前面所述的 L-BFGS 學習法則。另外，本研究所採用的轉換函數如下式所示：

⎪ ⎩ ⎪ ⎨ ⎧ − < − ≤ ≤ − > = 1 1 1 1 1 1 ) ( y if y if y y if y g (2.19) 值得一提的是，於訓練網路前，通常會將訓練數據加以正規化。因此，在整個網路運算過 程中，| y 值很少會大於 1。 | 確立類神經網路架構為本模式之重點所在。首先，吾人必須收集地震反應資料作為訓練 範例，加以各種合理之網路架構測試。在這過程之中，吾人需決定的重要因素有幾個地方： 1. 時間延滯的採用：所謂「時間延滯」，即式(2.17)與(2.18)中採用之m與n值。其在類 神經網路收斂與系統識別之準確度上有極重要之關係。一般而言，m與n若取的越大， 對本模式效果會越佳；但是相對的，網路訓練之時間將消耗越久，識別所得之準確度 提昇亦不大。 2. 隱藏層結點數的選擇：節點數之選擇依據問題而定，一般以試誤法來處理。在誤差容 許的範圍下，盡量減少其數目，以降低網路學習之時間。 3. 訓練資料數目的多寡：對類神經網路而言，資料之多寡即訓練案例數目之增減。一旦 有著較充足之資訊作分析，類神經網路可以在搜尋的空間中穩定地找到較佳解。然而 數目的多寡同樣主宰了網路學習的時間。一般情況下，常取資料中較具代表性之一段 資料作為訓練案例即可，如此既可不造成資料失真，亦可減少訓練的時間。 於是，根據類神經網路架構的程序，利用量測資料加以訓練之，以獲得適當的網路架構與 節點連接權值。而系統的動態特徵便可由這些連接權值以系統識別程序求得[38][39]。

2.4、模態參數估算程序

當完成網路訓練之後，網路的輸出輸入關係便可描述該系統的動力行為。依據本研究 之網路架構見圖(2.2），網路輸出與輸入之關係可近似地表示成如下：

{ }

([ ]{ } { }) ] ][ [ } {Y = W V X + W

θ

v +

θ

w _(2.20) 其中

{ }

T k t x t x t x

Y =(&&₁( ), &&₂( ),_L, && ( ))

；

{ }

T X =(XF) _(2.21)(2.22) 又 )), ( ) ( ) ( ) 2 ( ) 2 ( ) 2 ( ) 1 ( ) 1 ( ) 1 ( ( 2 1 2 1 2 1 m t x m t x m t x t x t x t x t x t x t x k k k − − − − − − − − − = && L && && L && L && && && L && && X (2.23) )) ( ) ( ) ( ) 1 ( ) 1 ( ) 1 ( ) ( ) ( ) ( ( 2 1 2 1 2 1 n t f n t f n t f t f t f t f t f t f t f l l l − − − − − − = L L L L F (2.24)

而矩陣

[ ]

W 與

[ ]

V 內之元素為wij與vjk_。向量

{ }

θw 與

{ }

θv 內之元素為θwi與θvj_{。將式(2.20)} 加以展開，可得

{ }

) ( ) ( ) ( ˆ ) ( ) ( ) ( ˆ ) ( ) ( ) ( 2 1 0 ) ( 2 2 1 1 ) ( 1 2 1 C j t f j t f j t f i t x i t x i t x t x t x t x l n j j k m i i k + ⎪ ⎪ ⎭ ⎪ ⎪ ⎬ ⎫ ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ − − − + ⎪ ⎪ ⎭ ⎪ ⎪ ⎬ ⎫ ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ − − − = ⎪ ⎪ ⎭ ⎪ ⎪ ⎬ ⎫ ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧

∑

∑

= = M && M && && && M && && W W (2.25) 其中

[ ][ ]

W V = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ 2 1 W Wˆ ˆ (2.26)

{ }

C =

[ ]

W

{ } { }

θ

v +

θ

w _(2.27) ] ˆ ˆ ˆ [ ˆ ( ) 1 ) 2 ( 1 ) 1 ( 1 1 n W W W W = _{L (2.28)} ] ˆ ˆ ˆ [ ˆ ( ) 2 ) 1 ( 2 ) 0 ( 2 2 m W W W W = _{L (2.29)} 式(2.25)與時間序列模式之 ARX 相似，而 ARX 可對等於結構系統的運動方程式。因此，結 構之動力特徵便可由 AR 之係數矩陣而求得[38]。根據網路之權值矩陣建立以下矩陣： ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = − (1) 1 ) 2 ( 1 ) 1 ( 1 ) ( 1 ˆ ˆ ˆ ˆ 0 0 0 0 0 0 0 0 0 0 0 0 ] [ W W W W I I I L M M M M M n n G (2.30) 於是從

[ ]

G 矩陣之特徵值與特徵向量，便可計算出系統之模態參數[39]。今假設λk與

{ }

ψk 分 別表示

[ ]

G 的第k個特徵值與特徵向量。而λk通常為複數，可表示為a_k+ib_k。則相對應的自 然頻率與振態阻尼比可經由下列式子計算之。 2 2 ~ k k k

σ

γ

γ

= + (2.31) k k k

σ

γ

ξ

=− /~ _(2.32) 其中γ~k_{為擬自然振動頻率，而}ξk_{為振態阻尼比。}γk與σk則根據下列式子決定。 ) ( tan 1 1 k k k a b t − Δ =

γ

(2.33)

) ln( 2 1 2 2 k k k a b t + Δ =

σ

(2.34) 其中1Δt_{為測量時的採樣頻率。} 由於

[ ]

G 矩陣本身之特殊結構，使得其特徵向量之間有以下的關係：

{ } { }

{ }

{ }

{ }

T T k n T k T k k T k k ( , , k ,..., k 1) 1 1 2 1 1 Ψ Ψ Ψ Ψ = Ψ

_λ

_λ

_λ

− (2.35) 其中

{ }

Ψk 1為相對於自然頻率γ~k之複數模態。 以上的推導過程中，乃基於類神經網路僅一層的隱藏層。然而，吾人只要將式(2.25)中 ) ( 1 ˆ i W 、 ) ( 2 ˆ j W 、及{C 稍作修正亦可將上列公式推展至多層隱藏層的情況。 }

2.5、類神經網路之預測誤差

當結構系統在外力作用下仍為線性，且其模態參數無明顯改變時，可根據前幾段時刻 所測量的反應及輸入，建立一個符合線性結構系統的類神經網路識別架構，如此將能夠準 確地，預測目前的結構反應。然而，若此時結構系統有所破壞或是退化的情形發生，則結 構系統便會呈現出非線性的行為。若以一健康、線性結構所建立的網路來預測此時的結構 反應，便會導致有較大的預測誤差。由 Masri 等人[13]所提出的指標可作為量化這些誤差 的工具。而為了簡單呈現，本報告僅採用各自由度預測值與實際量測值之間的平均絕對誤

差（mean absolute error, MAE），作為網路預測誤差的表示。MAE之定義如下：

MAE(i)= 1

1

( )

( )

T im ip t

y

t

y t

T

=

−

∑

(2.35) 其中i 表示第 i 個自由度，而yim與 yip分別為第i 個自由度上的正規化量測反應與預測反 應。

參、希伯特黃轉換於系統識別之基礎理論與架構

近年來由黃鍔博士等人發展出一種新的的訊號處理方法，它可以處理頻率域與時間域 的資料。該方法首先是利用經驗模態分解法(EMD)將訊號分解成數個內建模態函數，如此一 來每個內建模態函數將會是良好的希伯特轉換對，接著我們對每個內建模態函數進行希伯 特轉換，如此將可以得到頻率域與時間域的分解訊號，這樣的過程我們稱之他為希伯特黃 轉換。該方法是可適用於非穩態的訊號，而吾人將嘗試以此方法為基礎來分析多自由度的 結構中的自然頻率以及阻尼比。在利用希伯特黃轉換於多自由度線性系統的識別的過程 中，自由振動(Free Vibration)的訊號是吾人所需要的，我們第一可以利用衝擊載重可以 直接求得結構物自由振動的反應，第二若於現地的風力或地震力外擾的情況下，吾人必須 搭配使用隨機遞減法來求得各自由度的自由振動訊號。接下來可以利用帶通濾波器的技術 來處理得到的訊號，如此將可分離求得每個模態反應，將各個模態的反應透過經驗模態分 解法，分別求出每個模態的第一個內建模態函數，將此內建模態函數經由希伯特轉換後可 以得到其瞬時相位角與瞬時震幅對時間的關係曲線，最後使用最小平方法來近似求得的瞬 時相位角與瞬時震幅的對數值，透過簡單的計算，就可以求得多自由度系統中每個模態的 自然頻率與阻尼比。本章前幾節將介紹希伯特黃轉換相關理論，與隨機遞減法，之後為運 用在系統識別的架構。

3.1、簡介希伯特黃轉換

希伯特黃轉換(HHT)為1998年美國太空總署(National Aeronautics and Space Administration, NASA)的黃鍔博士所發表，其方法第一步是先將實驗所得之歷時訊號由經 驗模態分解法，過程是以不同的內部資料時間為尺度將能量直接萃取，如此訊號分解而產 生數個內建模態函數，而被分解出來的訊號就可以為非線性或非平穩性，如此一來這些分 量依就保留了原來訊號的物理特性，因此該分析方法對非線性(Non-linear)或非穩態 (Non-stationary)之歷時訊號亦有不錯的解析。經過分解後的每個內建模態函數將是良好 的希伯特轉換對，吾人也將可由此得到其瞬時頻率及振幅隨時間之分佈，使訊號能表現瞬 時變化的訊息與特性，稱為希伯特黃轉換。在黃鍔博士所提出的經驗模態分解法曾驗證過 該法具有完整性、正交性、局部性及可適性等之特性，因此也使得其在訊號處理中有著極 大的優勢。圖(3.1)所示為希伯特黃轉架構之流程圖。 1. 完整性：確保展開之精確度。 2. 正交性：保持能量的正定及避免能量不保守。 3. 局部性：使訊號表現瞬時變化的特性。 4. 可適性：適應暫態、非線性及非穩定性之系統條件。

3.2、即時頻率(Instantaneous Frequency)與希伯特轉換

在介紹經驗模態分解法前，須先瞭解何謂內建模態函數，而在說明內建模態函數前必

須先明確的定義即時頻率。即時頻率為分量函數所必要滿足的限制條件，而符合這些限制 條件的函數，就是所謂的內建模態函數。經驗模態分解法的目的即將原始訊號分解成多個 內建模態函數。 即時頻率的意義一直以來是廣受討論，1995年Cohen[40]提出了一個特殊單一分量訊號 的觀念，才使即時頻率有了一個具體的概念。在傳統的傅立葉分析中對於即時頻率的定義 以正弦及餘弦調和函數來做為基底，並且振幅是固定不變，即時頻率表示每一瞬間的頻率 值。就傅立葉分析而言，因為取樣的訊號必須大於一個完整的震盪周期，因此最少要有一 個完整的正弦波或餘弦波震盪周期來定義局部的頻率值。如此一來，對於頻率值隨著時間 改變的非穩態訊號而言是不合理的，因此我們藉助希伯特轉換將使訊號得以解析。對任意 的時間序列

X t

( )

，可以得到它的希伯特轉換

Y t

( )

，如下： 1 ( ') ( )= ' ' X t Y t p dt t t π ∞ −∞ −

∫

(3.1)

這裡的P 是指柯西主值(Cauchy Principal Value)，這個轉換對所有LP

階的函數皆存在 [41]。由此一定義，X(t)和Y(t)可以組合成一共軛複數，如此我們可以得到一解析的的訊號 Z(t): ) ( 2 1 2 2_{( ) ( )] ( )} [ ) ( i t e t a t Y t X t z = + = θ (3.2) 其中 1 2 2 ₂ ( ) [ ( ) ( )] a t = X t +Y t (3.3) ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ = − ) ( ) ( tan ) ( 1 t X t Y t θ (3.4) 理論上來說，有無限多種方法可來定義虛部，但是希伯特轉換提供了一個唯一的方法來定 義虛部，所以結果是一解析的函數。簡單來說，在式(3.1)中，希伯特轉換即定義為X(t)與 1/t的迴旋積分(Convolution)，因此，希伯特轉換的物理意義在於強調X(t)的局部特性。而 在式(3.2)中，利用極座標的表示式，更進一步的闡明了它局部性的本質，它將振幅及相位 角改變的三角調和函數利用最佳局部近似方法來近似X(t)。但即使有了希伯特轉換，對於式 (3.5)所定義即時頻率仍然是備受爭議的。

( )

d

t

dt

θ

ω

=

(3.5) 這促使了Cohen(1995)[40]提出“單一分量函數＂(Monocomponent Function) 這個專門的 名詞。在這理論中，對於資料將有一些必要的限制。在任何時間點上式(3.5)定義的即時頻 率必要是單一的值，它將只能表示一個分量，稱之“單一分量函數＂。不幸地，對於如何 判別是否為“單一分量函數＂的方法卻沒有提供一個清楚的定義。由於缺乏一精確的定 義，在文獻裡是採用“窄頻寬＂(Narrow Band) 來定義即時頻率。而為了得到有意義的即 時頻率，將使用Gabor[42]與Bedrosian[43]所討論而得的嚴格條件：對於任一個函數要得 到其有意義的即時頻率，它的傅立葉轉換的實部必須只能有正的頻率。對於資料分析而言，

要使這個條件成為具體可以實行的步驟因而發展出一個簡單的方法來加以應用。為了達成 這個目的，必須修正這個限制條件從全面性的變成局部性的，而且它的基底還必須滿足前 一節所列的幾個必要條件。

3.3、內建模態函數

前一節定義了即時頻率，其給予分量限制條件，為了使所有的分量都能滿足限制條件， 接下來就使用這些限制條件來定義內建模態函數。在物理上要定義一個有意義的即時頻率 其必要條件是函數的相對局部零均值(local zero mean)是否對稱的，並且跨零

(zero-crossings)的數量必須相等於極值(Extrema)數量總和。根據這些條件，定義了內建 模態函數。一個內建模態函數它是一個必須滿足下列條件的函數，其條件如下：

(1)在整個函數中，極值的數量必須要跟跨零點的數量相等或者相差一個；

(2)在任何時間點上，由極大值包絡線(The maxima value of envelope)及極小值包絡線(The minima value of envelope)所定義的均值包絡線(The mean value of envelope)必須是零。 其中這裡的極大值包絡線是由局部極大值所定義的，極小值包絡線是由局部極小值所定義 的。

第一個條件，與傳統平穩高斯過程(Stationary Gaussian Process)中窄頻寬的要求很 類似。第二個條件是一個新的建議，目的在將整體性要求改變成局部性的要求，這樣的好 處在於對處理不對稱的波形時即時頻率也將不會產生不必要的晃動；理想中，要處理訊號 的局部均值(local mean)應該要是零，但是對於非平穩性的訊號而言，局部均值是需要一 個局部的時間尺度(local time scale)來計算，這是非常不容易的。因此為了完成這個動 作，將利用局部極大值定義的極大值包絡線及由局部極小值定義的極小值包絡線來強迫局 部的對稱，來代替原來的方法，這樣的近似方法則可避免定義局部平均時間尺度的困擾。 內建模態函數它表達了深藏在訊號內部振盪模態，在每一個循環中由跨零點所定義的內建 模態函數只包含一個模態的振動，不會有很複雜的載波亦不會被限制到是一個窄頻寬的訊 號，它的振幅及頻率是可以變動的，頻率或振幅完全變動的函數皆可以為內建模態函數， 因此內建模態函數將可以是非平穩性的。綜觀之前的論述，若要得到即時頻率，需要解析 一個任意的資料成為內建模態函數分量，即時頻率才可以應用到每一個內建模態函數分量 裡去，對於複雜的資料來說，在同一時刻可能會得到超過一個以上的即時頻率。

3.4、經驗模態分解法

內建模態函數為希伯特黃轉換的前置條件，但大部分的輸入訊號都不是內建模態函數 這種形式，因為在任何給定的時間，可能包含了不只一個振動模態，這也就是希伯特轉換 不能對完整的訊號提供全盤頻率內涵的原因。因此必須分解輸入訊號成為多個內建模態函 數分量，要達到這樣的目的，則必須藉由經驗模態分解法，簡單的來說，經驗模態分解法 就是將輸入訊號分解成多個帶有物理特性的內建模態函數，而內建模態函數就如同傅立葉 轉換中的正弦、餘弦函數為希伯特黃轉換的基底。在經驗模態分解法中較特別的是其基底

函數是藉由分解原輸入訊號而得，這個方法是直觀的、直接的、後驗的、及可適性的。由 於分解的基底是從原來的訊號推導而來的，因此有下列假設條件: (1)欲分析之訊號至少有兩個極值，其中包括一個極大值和一個極小值; (2)訊號特徵時間尺度是定義成兩個極值之間的時間差; (3) 如果資料全無極值但只包含反曲點，可以將訊號微分一或多次，就 能將極值找出來。最後的結果可以由分量的積分得到。 根據經驗利用訊號中特徵時間尺度來定義其振動模態，然後依據它來分解訊號，這是一個 有系統的方法用來解析出內建模態函數，又可稱為轉移過程(Shifting Process)，此分解 法分別找出局部極大值和局部極小值，再將所有的局部極值利用立方雲線(Cubic Spline) 把它們連接起來，產生上圍的極大值包絡線和下圍的極小值包絡線(mean envelope)。則待 分析之原始訊號會被上下圍包絡線完全包覆，再由極大值包絡線與極小值包絡線取平均得 均值包絡線為 m1，而原始訊號與均值包絡線之差即是第一個分量，稱為h1，即： 1

( ) -

1

h

=

x t

m

(3.6) 其中x(t)是原始訊號，過程見圖(3.2)～(3.14)。此時的h1，如圖(3.14)，不一定是內建模 態函數，因原始資料突高及突低，會在如此轉移的過程中產生出新的極值或是平移誇大先 前存在的極值。而因為是直接利用原始資料取其極值包絡線，所以對極值包絡線的影響是 直接的，而對均值包絡線的影響是間接的，即使立方雲線近似再完美，斜率輕微的突起也 會放大成一個局部極值，使得局部零值從一個直角的座標系統改變成曲線的座標系統，換 句話說已經把均值包絡線和y軸當成新的座標系統。同時這個方法產生的新極值會使試驗中 消失的模態顯露，而靠著一連串的轉移也將會使得低振幅的載波復原。經過一次轉移過程 出來的結果將不能保證所有的局部極大皆為正，局部極小值皆為負，將不符合我們所需要 的內建模態。內建模態函數的導出跟原始訊號有很大的關係，所以當原訊號越複雜時，則 要達到符合內建模態函數條件的函數就需要較多次這樣的轉移過程。 對非線性的資料而言，EMD 仍然存在其他複雜的問題，如包絡線的均值可能會跟真正 的局部均值不一樣，所以，無論對訊號轉移多少次，一些不對稱的波形仍然會存在。除了 這些理論上的問題，在實際數值分析中，立方弧線近似在邊界的近似也會產生很嚴重的問 題，因為立方弧線近似會造成相當大的擺盪。更不理想的情形是末端的擺盪最後會傳遞到 訊號內部而將整個訊號破壞，我們可以試著在原始資料的前後各延伸一個餘弦波來消除末 端效應，便可求得改良的立方弧線近似，但也有其它的方法可以選擇。即使有這些問題， 轉移過程仍是可從資料中解析出必要的分量。 在第二次轉移過程中，則是將h1當作原始訊號，找出h1的極大值、極小值包絡線，然後 求得 h1的均值包絡線m11 ，即 11

1

-

11

h

=

h

m

(3.6) 如圖(3.15)～(3.18)，圖(3.18)為經過兩次轉移過程的結果，而經過兩次轉移過程所 得到的h11 將更加對稱，如果仍然有局部極大值並非為正，將必須要再繼續做轉移過程，圖 (3.19)即經過k次轉移過程的結果h1k ，已經完全符合一個內建模態函數的條件，以數學式

可表示成 1k

-

1(k 1) 1k

h

=

h

₋

m

(3.7) 圖(3.21)即是原始的訊號中第一個內建模態函數分量，利用轉移過程有兩個目的：(1) 消除載波；(2)使得波形振幅更為平順。藉著轉移過程消除載波以及平滑不平整的振幅，如 此過程才能產生有意義的即時頻率。振幅平滑是必需的，但將轉移過程次數過多，雖然平 滑了不平整振幅，但卻將訊號原先的物理特性消除，且通過轉移的次數達到一個極值，將 使得整個訊號變成固定振幅，而單純頻率變化的訊號；因此必須決定轉移過程停下來的收 斂條件。歷來，黃鍔博士等人曾提出兩種不同的收斂停止條件： (1) 前後轉移函數間的標準差小於我們預先設定的值，標準偏差(standard deviation, SD)公式定義如下：

( )

( )

2 1 0 2 1 0 T k k t k T k t

h

t

h t

SD

h

− = − =

−

=

∑

∑

(3.8) (2) 選定一個轉移次數值S，使得跨零點數量和局部極值總數量(局部極大值數量加上局 部極小值數量)要相等或差一個。 當滿足以上任一個收斂條件，即停止此一模態的轉移程序，再繼續下一個模態的轉移。在 第一個停止條件中而SD的大小該以多少為佳呢？S值取多少適宜呢？在黃鍔博士先前發表 的文章[1]中提出一般在SD值取0.2～0.3之間就已足夠，而S值的大小的範圍，也約介於4 ～8間就足夠了。 假定選擇了其中的一個停止條件，第一個內建模態將可因此產生。整體來說， C1應該 是包含訊號中最佳的時間尺度或者是最短週期的分量。從原來訊號中分離C1可得到：

( )

1 1

r

=

x t

−

c

(3.9) 因為餘數r1，仍然包含較長週期變動的分量，所以它即被當成新的資料並且再利用上述相 同的轉移過程來做處理。這個過程持續的重複在之後所有餘數rj，公式如下： 2 1 2 1 n n n

r

r

c

r

r

₋

c

= −

=

−

M

(3.10) 轉移過程可以靠著任何下面所提的準則來停止，一是當分量Cn或者是rn變得比預定值來的 小時就停止，二是當餘數rn，變成單調的函數，使得沒有IMF 可以解析出即停止。即使是 零均值的訊號，最後的餘數仍可能不是零;假如該訊號有一個趨勢(trend)，轉移過程後最 後產生的餘數就是其趨勢。將Eq(3.9)與(3.10)加起來，最後可以得到

( )

1 n j n j

x t

c

r

=

=

∑

+

(3.11)

如此一來，將可以把一筆資料分解成n個經驗模態(empirical modes)及一個可以當作均值

趨勢(mean trend)或常數的餘數rn。當應用EMD時，不需要任何的均值或零值參考軸，只需

要知道局部極值的位置。在轉移過程中，每個分量的零值參考軸即會自動產生。且由於EMD 不需要零值參考軸，故在對具有非零均值的資料時，EMD 即可免除資料中尋找均值的這項 麻煩的步驟。

3.5、隨機遞減法

假設一自由度為n之線性系統受到一零平均值的穩態白訊隨機振動(zero-mean

stationary white-noise random vibration)外力，其運動方程式可寫為：

( )

( )

( )

( )

.. .

M x t

+

C x t

+

Kx t

=

F t

　 (3.12) 其中M、C及K分別代表n個自由度的結構之質量、阻尼以及勁度矩陣，x(t)為位移反應，而 F(t)為外力向量。考慮此系統於某時刻ti滿足運動方程式，接著將時間平移t後可得：

(

)

(

)

(

)

(

)

.. . i i i i

M x t

+ +

t

C x t

+ +

t

Kx t

+ =

t

F t

+

t

(3.13) 將(3.13)式進而累加N 次後可導出：

(

)

(

)

(

)

(

)

.. . 1 1 1 1 N N N N i i i i i i i i

M

x t

t

C

x t

t

K

x t

t

F t

t

= = = =

⎡

₊

⎤

₊

⎡

₊

⎤

₊

⎡

₊

⎤

₌

₊

⎢

⎥

⎢

⎥

⎢

⎥

⎣

∑

⎦

⎣

∑

⎦

⎣

∑

⎦

∑

(3.14) 接著將(3.14)式除以N 則可寫成：

( )

( )

( )

(

)

.. . 1

1

N i i

M y t

C y t

Ky t

F t

t

N

=

+

+

=

∑

+

(3.15) 其中

( )

(

)

.. .. 1

1

N i i

y t

x t

t

N

=

=

∑

+

；

( )

(

)

. . 1

1

N i i

y t

x t

t

N

=

=

∑

+

；

( )

(

)

1

1

N i i

y t

x t

t

N

=

=

∑

+

(3.16) 由於假設F(t)為具零平均值之穩態白訊，因此當N→∞時，(3.15)式等號右邊必為零向量， 而(3.15)式可改寫成：

( )

( )

( )

.. .

0

M y t

+

C y t

+

Ky t

=

(3.17) 如此即可把原本(3.12)式之外力運動方程式轉化為(3.17)式的自由振動方程式，只是變數 向量由x(t)改成y(t)。一般而言，當(3.16)式之訊號疊加超過500次(亦即N>500)時，即可

得穩定的特徵曲線y(t) [44]，文獻中稱為隨機遞減曲線(random decrement signature)。

如圖(3.22)所示之任意速度訊號為例，隨機遞減曲線之求取可分為以下步驟詳細說明：

為基準，使得x t.

( )

在時間為t1、t2一直到tN 等N 個時間點通過 .. sd x 值，如圖(3.22)所示。 如此在這些時間點， x t.

( )

之斜率將依序正負相間。x 取訊號平均值(零)加上標準差..sd (standard deviation) ，是為了確保疊加次數足夠使其疊加後之反應接近自由振動。 2. 訂定截取訊號的延時Td ，通常Td 須取為系統預估最大自然週期的數倍[45]，如此可以 使最後萃取的訊號能明顯地表示結構之動態特性。 3. 由t1 、t2等開始依序截取Td 時間長之量測訊號，相繼累加N 次而產生隨機遞減曲線：

( )

(

)

.. .. 1

1

, 0

N i d i

y t

x t

t

t

T

N

=

=

∑

+

< <

(3.18) 由於

( )

.. .. sd i

x t

=

x

，i=1,2,3,…,N (3.19) 所以

( )

.. ..

0

sd

y

=

x

(3.20) 另外

( )

...

0

i

x t

≥

，i=1,3,5,7…,2n-1 (3.21) 且

( )

...

0

i

x t

≤

，i=2,4,6,8,…,2n (3.22) 上述步驟隨著疊加次數的增加，因初始加速度所引起的自由振動反應將由於相鄰兩段疊加 截取訊號之起始斜率大小相似、正負相反而互相抵銷。另外，由(3.20)式可以清楚知道隨 機遞減曲線y(t)乃代表僅由初始加速度 .. sd

x

所引起的自由振動反應。

3.6、多自由度系統之模態反應

如同前節，一個n自由度為之結構其運動方程式可寫為：

( )

( )

( )

( )

.. .

M X t

+

C X t

+

KX t

=

F t

 (3.23) 其中

X t

( )

=

[

x x

₁

, ,...,

₂

x

_n

]

T=n 個自由度的位移向量，F t

( )

=n 個自由度的外力向量，Ｍ，Ｋ， Ｃ為(n×n)的質量、勁度、阻尼矩陣。假設有振動模態產生，位移加速度反應可以分離成 n 個真實的模態：

( )

( )

1 n j j j

X t

q t

=

=

∑

Φ

；

( )

( )

.. .. 1 n j j j

X t

q t

=

=

∑

Φ

(3.24) 上式中，

Φ

_j是第 j 個模態向量，表是第 j 個振態的形狀(mode shape)，而q t_j

( )

是第 j 個 模態的廣義模態座標(generalized modal coordinate)。將(3.24)代入(3.23)並依據振態 形狀的正交性，我們可以把(3.23)分離成 n 個模態的動力方程。

( )

.. . 2

2

T

/

j j j j j j j j

q

+

ζ ω

q

+

ω

q

= Φ

F t

m

(3.25) 其中

ω

_j為第 j 個模態頻率，

ζ

_j為第 j 個阻尼比，

m

_j為第 j 個模態的質量。 假設若有個衝擊載重施加在第 k 個自由度上，如：

f t

_k

( )

=

F

₀

δ

( )

t

，其他 k 之外的自 由度皆為零

f t

_j

( ) 0

=

，

f t

_j

( )

是表示外力向量

F t

( )

矩陣中的第 j 個元素的值。然後，第 j 個廣義模態座標的加速度反應可以得到如下： .. 0 2

( )

cos

2

1

j jt kj j dj j j j j

F

q t

e

t

m

ζ ω

φ ω

_ω

_ϕ

π

ζ

−

⎛

⎞

=

_⎜

−

+

_⎟

⎝

⎠

−

(3.26) 其中

φ

_kj表示第 j 個振態向量Φ 矩陣中第 k 個元素，_j

ω

_dj為第 j 個模態的阻尼頻率， 1 2 2 2

tan

2

1

/(1 2 )

j j j j

ϕ

₌

−

⎡

ζ

₋

ζ

₋

ζ

⎤

⎣

⎦

為第 j 個模態的相位差。而在第 p(p=1,…,n)個自由 度的加速度反應 ..

( )

p

x t

可表示如下： .. .. .. 1 1

( )

n

( )

n

( )

p _{pj j} pj j j

x t

φ

q t

x

t

= =

=

∑

=

∑

(3.27) 其中 .. .. , ,

( )

( )

cos

2

j jt pj _{pj j pj k dj j pj k}

x

t

=

φ

q t

=

B

e

−ζ ω

_⎜

⎛

ω

t

−

ϕ

+ +

π

ϕ

⎞

_⎟

⎝

⎠

(3.28) 0 , ₂

1

pj kj j pj k j j

F

B

m

φ φ ω

ζ

=

−

(3.29) 在(3.28)，

ϕ

_{pj k,} 表示在第 j 個振態中，第 p 個元素與第 k 個元素的相位差。假若存在簡正 模(normal mode)，則所有的振態形狀都將會是實數，如此

ϕ

_{pj k,} 將會是

±

2m

π

或者為

2(

m

1)

π

±

+

，m 為整數，如： 當

ϕ

_{pj k,}

= ±

2

m

π

，則

φ φ

_pj

/

_kj

>

0

當

ϕ

_{pj k,}

= ±

2(

m

+

1)

π

，則

φ φ

_pj

/

_kj

<

0

若把雜訊的干擾也列入考量的話。則量測出來的加速度反應向量 .. .. .. .. 1 2

( ) [ ( ), ( ), , ( )]

T n

Z t

=

z t z t

_L

z t

將可表示如下： .. ..

( )

( )

( )

Z t

=

X t

+

V t

(3.30) 其中

V t

( ) [ ( ), ( ), , ( )]

=

v t v t

_{1 2}

K

v t

_n T為一個噪訊向量。量測出來第 p 個自由度受衝擊載重的 加速度反應 ..

( )

p

z t

將可表示如下： .. .. .. 1

( )

( )

( )

n

( )

( )

p p _p pj _p j

z t

x t

v t

x

t

v t

=

=

+

=

∑

+

(3.31) 其中 ..

( )

pj

x

t

在(3.26)到(3.29)式中已清楚說明，

v t

_p

( )

我們將試著以零均值高斯白噪作為模 擬。

3.7、模態參數之識別程序

經驗模態分解法可以分解量測的反應訊號 ..

( )

p

z t

成為適於希伯特轉換的內建模態。依 據經驗模態分解法且搭配著帶通濾波器，可以從 ..

( )

p

z t

中分解出每個模態的反應。應用帶通 濾波器的有點有兩個優點： (1) 便於處理訊號擁有較高模態頻率的訊號， (2) 去除訊號受到高頻的雜訊干擾。 把加速度反應 ..

( )

p

z t

經過傅立葉轉換得到的傅立葉頻譜圖中，可以看到每個自然頻率大 致的頻帶分佈情形，假設得到第 j 的模態的訊號頻率

ω

_j介於

ω

_jL與

ω

_jH之間如：

(

1,2, , )

jL j jH

j

n

ω

<

ω

<

ω

=

_K

，接著將量測到的訊號 ..

( )

p

z t

通過

ω

_jL

<

ω

_j

<

ω

_jH的帶通濾 波器，將可分別得到這個模態的歷時訊號，接著再利用經驗模態分解法處理之，經過大量 的轉移次數，所得到的第一個內建模態，將相當接近第 j 個模態的反應 .. pj

x

。重複處理

1,2, ,

j

=

_K

n

將可以得到這 n 個模態的反應 ..

(

1,2, , )

pj

x

j

=

_K

n

。這個方法可以很輕易的從 ..

( )

p

z t

分離出每個模態的反應，且可以移除這個頻段

(

ω

_jL

<

ω

_j

<

ω

_jH

)

以外的雜訊。不過在 應用帶通濾波器時，盡可能的讓相位偏移(phase shift)越小越好。 經由上面的過程可以得到模態反應 .. pj

x

，接著使用希伯特轉換，可得到 ~

( )

pj

x

t

&&

如下所示： .. ~ .. _{( )} ( ) [ ( )] ( ) pj pj pj x t x t HT x t d t

τ

π

τ

∞ −∞ = = −

∫

&& (3.32) 而由(3.28)知道 .. pj

x

： .. .. , ,

( )

( )

cos

2

j jt pj _{pj j pj k dj j pj k}

x

t

=

φ

q t

=

B

e

−ζ ω

_⎜

⎛

ω

t

−

ϕ

+ +

π

ϕ

⎞

_⎟

⎝

⎠

應用 Bedrosian 定理將能得到 ~

( )

pj

x

t

&&

如下： , , , , , ~ ( ) ( )sin ( )cos 2 2 pj _{pj k LP j dj j pj k HP j dj j pj k} x t =B _⎢⎡a t _⎜⎛ω t−ϕ + +π ϕ ⎞_⎟+a t ⎛_⎜ω t−ϕ + +π ϕ ⎞_⎟⎤_⎥ ⎝ ⎠ ⎝ ⎠ ⎣ ⎦ && % (3.33) 由(3.29)式知其中 _, 0 2

1

pj kj j pj k j j

F

B

m

φ φ ω

ζ

=

−

，而 , ₀ 2 2 2

2

1

( )

dj _{j j}

cos( )

LP j j j

a

t

ω

ζ ω

ω ω

t d

π

ζ ω

ω

=

+

∫

； ~ , _{2 2 2}

2

1

sin( )

dj j j HP j j j

a

t d

ω

ζ ω

ω ω

π

ζ ω

ω

∞

=

+

∫

(3.34) 得到了 .. pj

x

與 ~

( )

pj

x

t

&&

，如此第 j 個模態的解析訊號將可表示如下： .. ~ ( )

( )

( )

( )

( )

i pj t pj pj pj pj

Y

t

=

x

t

+

i x

&&

t

=

A t e

θ (3.35) 其中的瞬時振幅

A

_pj與相位角

θ

_pj如下：

2 2 , , 1 ~ 2 2 , , , ,

( )

{

cos

2

[

( )sin

( )cos

] }

2

2

j jt pj pj k dj j pj k HP j LP j dj j pj k dj j pj k

A t

B

e

t

a

t

t

a

t

t

ζ ω

_ω

_ϕ

π

_ϕ

π

π

ω

ϕ

ϕ

ω

ϕ

ϕ

−

⎛

⎞

=

_⎜

−

+ +

_⎟

+

⎝

⎠

⎛

₋

_{+ +}

⎞

₊

⎛

₋

_{+ +}

⎞

⎜

⎟

⎜

⎟

⎝

⎠

⎝

⎠

(3.36) ~ 1 , , ,

( ) tan {

[

( ) tan

( )]}

2

j jt HP j pj

t

e

a

LP j

t

dj

t

j pj k

a

t

ζ ω

π

θ

₌

−

⎛

ω

₋

ϕ

_{+ +}

ϕ

⎞

₊

⎜

⎟

⎝

⎠

(3.37) 假設

ζ

_j遠小於

ω

_j的狀況下，式(3.34)將可改寫如下 , ₀ 2 2 2

2

1

( )

dj j j

cos( )

_{j j}t LP j j j

a

t

ω

ζ ω

ω ω

t d

e

ζ ω

π

ζ ω

ω

−

=

≈

+

∫

； ~ , _{2 2 2}

2

1

sin( )

0

dj j j HP j j j

a

t d

ω

ζ ω

ω ω

π

ζ ω

ω

∞

=

=

+

∫

(3.38) 因此(3.33)將變成 ~ , ,

( )

sin

2

j jt pj _{pj k dj j pj k}

x

t

=

B

e

−ζ ω

⎛

_⎜

ω

t

−

ϕ

+ +

π

ϕ

⎞

_⎟

⎝

⎠

&&

(3.39) 瞬時振幅

A

_pj與相位角

θ

_pj也將可改寫為 ,

( )

j jt pj pj k

A t

=

B

e

−ζ ω (3.40) ,

( )

2

pj

t

dj

t

j pj k

π

θ

=

ω

−

ϕ

+ +

ϕ

(3.41) 將(3.40)取自然對數後再對時間 t 微分則：

ln

_pj

( )

j j

A t

dt

= −

ζ ω

(3.42) 而將(3.41)對 t 微分則：

( )

( )

pj pj dj

d

t

t

dt

θ

ω

=

=

ω

(3.43) 1 2 2

(1

)

dj j j

ω

=

ω

−

ζ

(3.44) 因此在

ζ

_j遠小於

ω

_j的狀況下，可以利用上面三式計算出每個模態的自然頻率

ω

_j與阻尼比

j

ζ

。但若在

ζ

_j並非很小的狀況下，必須重新考慮(3.36)與(3.37)二式，其中

ln

A

_pj與

θ

_pj將 不在是對時間 t 的線性方程式，而黃鍔博士等人隨後做出的數值驗證結果[22]發現到，振 幅的變化會導致頻率會在均值附近產生波動，且不會改變均值，因而提出可以使用最小平 方法近似的方式來估算自然頻率與阻尼係數。