[發明專利]一種基于狀態觀測器的振動模態參數識別方法

有效

• 申請號: 201710548811.9
• 申請日: 2017-07-07
• 公開（公告）號: CN107391818B
• 公開（公告）日: 2019-10-11
• 發明（設計）人: 盛賢君;周少征
• 申請（專利權）人: 大連理工大學
• 主分類號: G06F17/50
• 分類號: G06F17/50;G01M7/02;G01H17/00
• 代理公司: 大連理工大學專利中心 21200
• 代理人: 關慧貞

1.一種基于狀態觀測器的振動模態參數識別方法，其特征是，該方法針對特定結構和關注的振動階次，先利用ANSYS模態分析近似得到各階振動頻率和振型，再建立自由振動實驗，采用多個非接觸式的高分辨率激光位移傳感器測量結構位移信號，據此建立整個系統的模態狀態空間方程和狀態觀測器；在利用狀態觀測器得到各階自由振動模態坐標之后，采用FIR濾波器濾除其中的干擾信號；最后采用指數衰減的正弦函數來擬合濾波后的模態坐標獲取各階振動精確的頻率和阻尼比，利用錘擊法在自由振動實驗中獲得更精確的模態參數；方法的具體步驟如下：

步驟一：組建自由振動實驗系統，

自由振動實驗系統由鋁懸臂梁(1)、鋁配重(2)、鋼配重(3)、固定端(4)、激光位移傳感器組(5)、數據采集卡(6)和計算機(7)組成；鋁懸臂梁(1)的左端為固定端(4)，鋁懸臂梁(1)的右端裝有鋁配重(2)和鋼配重(3)，鋁懸臂梁(1)的上表面裝有激光位移傳感器組(5)，激光位移傳感器組(5)、數據采集卡(6)與計算機(7)連接；以激光位移傳感器組(5)作為測量元件測量懸臂梁的位移信號，并經數據采集卡(6)進行A/D轉換后傳輸給計算機(7)，在計算機(7)內對離散化的位移信號進行數據處理得到精確的懸臂梁模態參數；

步驟二：建立結構的幾何模型，確定結構的材料力學參數，并導入到ANSYS中進行模態分析；

首先確定結構的材料力學參數，并導入到ANSYS中進行模態分析，分析結束后，將關注的振動階次的頻率(ω₁,ω₂,…,ω_n)和正則振型導出；選取自由振動實驗中激光位移傳感器的測量點和測量方向，激光位移傳感器的個數選為關注的振動階次的個數n；

將各階振動阻尼比取為0，系統的模態阻尼矩陣D和模態剛度矩陣K為：

系統的模態坐標動力學方程為：

選取n維振動模態坐標向量η和其變化率為狀態量，建立模態狀態空間方程如下：

其中，x為2n維狀態向量，y為各傳感器測量組成的n維輸出向量；A為系統矩陣，為2n階方陣；C為輸出矩陣，為n行2n列矩陣：

式中：I_n×n為n階單位陣；C_η為n階方陣，其第i行第j列元素其中，和分別代表第i個激光位移傳感器測量點在結構幾何建模坐標系中的位置和測量方向，代表結構的第j階正則振型矢量場；

在建立完狀態空間方程后對系統的能觀測性進行校驗，如不滿足，重新選取傳感器的測量點，直到滿足能觀測性要求；系統完全能觀測的充要條件是能觀測性矩陣N滿秩；N可以由系統矩陣A和輸出矩陣C構造出來：

步驟三：進一步以y和的組合作為反饋量建立狀態觀測器方程：

式中：

是x的狀態估計值，狀態觀測器的增益矩陣L由MATLAB中極點配置的函數來確定，具體程序為：

L＝(place(A',Cz',V))'； (8)

其中，V是由期望的狀態觀測器極點所構成的向量，極點選取為對應階次圓頻率的2～3倍的負數可以保證狀態觀測器快速跟蹤傳感器信號變化；

步驟四：利用錘擊法做自由振動實驗，采集各傳感器的測量數據，經過狀態觀測器得到各階振動模態坐標隨時間的變化；針對其中的雜波成分，設計FIR帶通濾波器將其濾除，FIR濾波器的離散時域表達式為：

其中，為濾波器輸入信號序列的第i-k個值，為濾波后的輸出信號序列的第i個值，N為的抽樣數，一般抽樣數越多，濾波器越接近理想濾波器；h(k)為濾波器的單位脈沖響應序列h的第k個值，h由MATLAB中的FIR濾波器設計函數得到，具體程序為：

h＝fir1(N,[f1,f2]*2/Fs,'DC-0')； (10)

其中，f1和f2分別為帶通濾波器的低截止頻率和高截止頻率，根據ANSYS分析結果和估計誤差范圍進行設定；Fs為采樣頻率，應取為最高關注頻率的20倍以上；'DC-0'表示選取帶通濾波器設計模式；在得到h之后，利用以下MATLAB程序查看所設計的FIR濾波器的幅頻特性和相頻特性：

freqz(h,1)； (11)

步驟五：在對各階觀測到的模態坐標進行濾波之后，用指數衰減的正弦變化函數對其進行擬合；首先查找信號中的極大值點，對應的時間組成的數組為time，極大值組成的數組為local_max，然后對這些極大值點做指數函數擬合，用公式(12)、(13)、(14)、(15)：

f＝(numel(time)-1)/(time(end)-time(1))； (12)

myfunc＝inline('b(1)*exp(b(2)*x)','b','x')； (13)

b＝nlinfit(time,local_max,myfunc,[0 0])； (14)

c＝-b(2)/(2*pi*f)； (15)

則精確的f和阻尼比c利用MATLAB程序計算出來。