基于離散傅里葉變換的頻響函數(shù)非參數(shù)辨識

基于離散傅里葉變換的頻響函數(shù)非參數(shù)辨識

頻響函數(shù)是振動結(jié)構中輸出信號和輸入信號的傅里葉變換比率。振動結(jié)構的頻響函數(shù)包括原振動結(jié)構物理參數(shù)的所有信息,如模態(tài)參數(shù)的確定、振動結(jié)構損失的程度。反之振動結(jié)構物理參數(shù)的信息也唯一地決定了頻響函數(shù)的變化。目前的研究側(cè)重直接從頻響函數(shù)中辨識振動結(jié)構的模態(tài)參數(shù),相對于模態(tài)參數(shù)而言,頻響函數(shù)更直接、誤差小、含有更豐富的原始數(shù)據(jù)信息。正是頻響函數(shù)的這些優(yōu)勢,可將頻響函數(shù)作為模態(tài)參數(shù)辨識和振動結(jié)構監(jiān)測損傷的主要依據(jù)。雖對頻響函數(shù)的研究較廣泛,但絕大部分的研究都是假設已知的頻響函數(shù)事先存在,或利用原始定義直接得到頻響函數(shù)。文獻為了克服對輸入-輸出觀測序列作離散傅里葉變化時所帶來的暫態(tài)和泄露譜影響,將通常的離散傅里葉變化推廣到同時考慮初始和終端狀態(tài)。采用頻響函數(shù)數(shù)據(jù)的頻域多輸入多輸出狀態(tài)空間模型,用頻響函數(shù)數(shù)據(jù)代替?zhèn)鹘y(tǒng)的狀態(tài)空間模型中的輸出,進而采用子空間辨識法求解狀態(tài)空間模型中的系統(tǒng)矩陣。但文獻中的子空間辨識法能使用的前提條件是在各個頻率點處的頻響函數(shù)估計值在統(tǒng)計概率意義下達到漸近有效性和一致性。文獻分析振動結(jié)構在環(huán)境激勵下的模態(tài)參數(shù)辨識問題,其使用的辨識方法仍是子空間辨識法,以繞開頻響函數(shù)的不準確而使得直接由不準確的頻響函數(shù)得到較大偏差的模態(tài)參數(shù)值。文獻直接利用頻響函數(shù)識別結(jié)構損傷,即頻響函數(shù)曲率法的原理類似于振型曲率法,但不需要測試振型,比振型曲率法識別效果更佳。以上代表在振動工程領域?qū)︻l響函數(shù)的研究還僅僅停留在頻響函數(shù)的應用價值取向上,還未對頻響函數(shù)的本質(zhì)做進一步的理論研究。但在系統(tǒng)辨識領域已對頻響函數(shù)做了一些基礎性的鋪墊工作。如文獻提出了多種不用于原始定義的頻響函數(shù)估計表達式,并分析了這些表達式各自的統(tǒng)計性質(zhì)。文獻從二階統(tǒng)計量-方差的角度分析為得到漸近無偏的頻響函數(shù)估計值,需要在對輸入-輸出觀測序列的離散傅里葉變換過程中引入加權濾波器,以達到平滑濾波噪聲的功能。文獻分析線性動力系統(tǒng)中頻響函數(shù)的頻域極大似然估計。文獻分析非線性系統(tǒng)的最優(yōu)線性無偏估計,即用一個頻響函數(shù)估計值來最優(yōu)地逼近于原非線性系統(tǒng)。本文深入研究頻響函數(shù)的辨識估計,為振動工程領域中頻響函數(shù)的直接使用提供準確表達式,將頻響函數(shù)估計方法應用于振動結(jié)構的模態(tài)參數(shù)辨識之中。頻響函數(shù)的估計通常有兩種方法:參數(shù)法和非參數(shù)法。頻響函數(shù)估計的參數(shù)法中,當增加了額外的未知參數(shù)勢必造成參數(shù)辨識成為一個奇異問題,不能得到關于未知參數(shù)的唯一解。因此本文研究頻響函數(shù)估計的非參數(shù)法。根據(jù)頻響函數(shù)的定義可知,首要一步是對輸入-輸出觀測序列做離散傅里葉變換,在此變換過程中同樣考慮初始和終端狀態(tài)帶來的暫態(tài)泄露項和觀測噪聲譜項對頻響函數(shù)估計的影響。在增加此兩項的基礎上,為得到頻響函數(shù)估計值,聯(lián)合頻響函數(shù)、初始-終端狀態(tài)和脈沖響應系數(shù)作為整體的待辨識未知參數(shù)矢量,將頻響函數(shù)估計問題轉(zhuǎn)化為求解一個線性最小二乘優(yōu)化問題。針對此線性最小二乘優(yōu)化問題的特殊形式,提出一種可分離的求解過程1結(jié)構模態(tài)參數(shù)辨識頻響函數(shù)主要應用于模態(tài)參數(shù)辨識中,在時域環(huán)境下進行模態(tài)參數(shù)辨識,通常采用微分或差分方程、傳遞函數(shù)矩陣等形式來描述系統(tǒng)結(jié)構的振動行為。而以上輸入-輸出模型只刻劃了結(jié)構的外部特性,并未深入內(nèi)部。而狀態(tài)空間模型則可以深入到結(jié)構的內(nèi)部情況。此處介紹如何將輸入輸出模型轉(zhuǎn)化為狀態(tài)空間方程模型,以及得到的頻響函數(shù)與狀態(tài)方程模型之間的關聯(lián)。在振動工程領域經(jīng)常遇到振動響應問題:對于非時變、粘性阻尼、線性自由度為N的系統(tǒng),可以根據(jù)Hamilton原理,推導出在外力u(t)作用下的運動方程為:Μ??y(t)+C1˙y(t)+Κy(t)=u(t)(1)式中:M,C,K∈RN×N分別為振動系統(tǒng)的質(zhì)量矩陣、粘性矩陣和剛度矩陣,y(t)是位移向量。式(1)給出系統(tǒng)的廣義運動方程,且該式是結(jié)構動力學研究最常用的形式,也是結(jié)構模態(tài)參數(shù)辨識的基本方程式。將上式與某些恒等式組合在一起可得:˙y(t)=Ι˙y(t)??y(t)=-Μ-1Κy(t)-Μ-1C1˙y(t)+Μ-1u(t)}(2)假設狀態(tài)變量x(t)為2N維矢量:x(t)=[y(t)˙y(t)]Τ可得狀態(tài)方程為:˙x(t)=[˙y(t)??y(t)]=A1x(t)+B1u(t)(3)其中:A1,B1分別為:A1=[0Ι-Μ-1Κ-Μ-1C1]?B1=0Μ-1Ι其中:A1∈R2N×2N為系統(tǒng)矩陣,B1∈R2N×1為輸入矩陣。系統(tǒng)的輸出向量與狀態(tài)向量x(t)間的關系為:y(t)=Cx(t)+v(t)(4)其中:y(t)為m維的輸出向量,C為m×2N維的結(jié)構輸出矩陣,v(t)為m維的系統(tǒng)測試輸出噪聲。聯(lián)合(3)、(4)式所述狀態(tài)方程的離散化形式為:xn+1=Axn+Bunyn=Cxn+vn}(5)通常對(5)式進行離散傅里葉變換可得,此時經(jīng)常忽略初始狀態(tài)和終端狀態(tài)下,系統(tǒng)的頻響函數(shù)矩陣為:G(z)=C(zI-A)-1B(6)求取狀態(tài)矩陣A的特征值和特征向量,則有如下等式成立:A=PΛP式中:Ρ=[φ(1)φ(2)?φ(2Ν)]Τ,且φ(i)為A的特征向量,Λ=diag(λ1λ2?λΝ),且λi為A的特征值。代入可將狀態(tài)空間模型轉(zhuǎn)化為模態(tài)展開模型:G(z)=CΡ(zΙ-Λ)-1Ρ-1B=[φφ*](zΙ-Λ)-1[LL*](7)上式中稱矩陣A的特征值為系統(tǒng)的極點,CP為模態(tài)的振型矩陣,φ-1P為模態(tài)的貢獻因子矩陣。以上的具體分析在于系統(tǒng)矩陣A的子空間辨識求解,但若能求得頻響函數(shù)G(z)的估計值?G(z),則可以利用頻響函數(shù)估計值?G(z)來辨識模態(tài)參數(shù)。頻響函數(shù)矩陣式(6)中第(i,j)個元素可表示在第j自由度方向施加單位簡諧荷載時,在第i自由度方向產(chǎn)生的位移響應。響應測點i和激勵點j之間的位移頻響函數(shù)為:Gij=Ν∑r=1φirφjrΚr-ω2rΜr+jωrCr(8)式中:ωr為振動結(jié)構的第r階模態(tài)頻率;Kr,Mr,Cr是對應的模態(tài)剛度、模態(tài)質(zhì)量和模態(tài)阻尼;φir,φjr分別是對應的(i,j)點處的振型幅值。由式(8)可知,若左邊的頻響函數(shù)估計值求出后,可將其改寫成分式和的形式,從各個分式中的系數(shù)即可直接得到對應各個模態(tài)參數(shù)值。2散傅里葉變換重寫式(5)所示的離散化狀態(tài)空間模型形式:xn+1=Axn+Bun,yn=Cxn+vn上式的頻響函數(shù)定義為:G0(ejω)=∞∑t=0g0te-jωt=∞∑t=1CAtBe-jωt(9)其中:g0t=CAtB表示系統(tǒng)的脈沖響應系數(shù),本文的研究目標在于,從一族輸入-輸出觀測序列對{u(n),y(n)}Νn=1(N表示觀測數(shù)據(jù)個數(shù))中,辨識出在所有頻率網(wǎng)格點ωk=2πΝk?(k=0?Ν-1)處的頻響函數(shù)G0(ejωk)的非參數(shù)估計值?G0(ejωk)。當N個輸入-輸出觀測數(shù)據(jù)對可獲得時,應用N點的離散傅里葉變化公式:X(k)=F(xn)=(1√Ν)Ν-1∑n=0xnz-nkX(k+1)=F(xn+1)=1√ΝΝ-1∑n=0xn+1z-nk=zkX(k)+zk√Ν(xΝ-x0)zk=exp(j2πk/Ν)}(10)式(10)關于xn+1的離散傅里葉變換過程中,考慮系統(tǒng)的初始狀態(tài)x0和終端狀態(tài)xN的存在,且推廣了式(6)中的一般廣義情況。將式(10)應用于式(5),并經(jīng)過整理可得:YΝ(ωk)=G0(ejωk)UΝ(ωk)+1√ΝΝ-1∑t=0ytratejωkt+VΝ(ωk)=G0(ejωk)UΝ(ωk)+ΤΝ(ejωk)+VΝ(ωk)(11)其中:ytrat=CAt(x0-xN)式中:x0為初始狀態(tài),xN為終端狀態(tài),第2項TN(ejωk)表示暫態(tài)項或泄漏譜,其產(chǎn)生的原因在于有限觀測時間下初始和終端狀態(tài)帶來的影響。且由TN(ejωk)的表達式可見,只有當用周期信號激勵系統(tǒng)時,有x0=xN,從而暫態(tài)項可消除。對于式(11)中的真實頻響函數(shù)G0(ejωk)的估計值求法通常采用如下定義等式:?G(ejωk)=YΝ(ωk)UΝ(ωk)+ΤΝ(ejωk)UΝ(ωk)+VΝ(ωk)UΝ(ωk)此估計值的特點在于,估計值?G(ejωk)僅定義在固定的頻率點處;此估計值是無偏估計,且其方差隨著1/N而遞減;因沒有平滑濾波的功能存在,該估計值的方差不會隨著N的增加而遞減;在隨機激勵作用下,該估計誤差將會任意大。因此根據(jù)輸出信號和輸入信號的離散傅里葉變化之比所得到的頻響函數(shù)估計是不理想的,不能直接應用于振動模態(tài)參數(shù)辨識之中。為了從輸入輸出觀測數(shù)據(jù)序列的離散傅里葉變換UN(ωk)、YN(ωk)中辨識在整個頻域范圍ωk=2πΝk?(k=0?Ν-1)中的頻響函數(shù)G0(ejωk)的估計值?G(ejωk)。3線性關系的建立因式(11)中不僅含有頻響函數(shù),還含有暫態(tài)項,同時將暫態(tài)項和頻響函數(shù)作為需要辨識的量。但式(11)展開后僅有N個復方程(k=0,…,N-1),而辨識過程中卻有3N個未知待辨識量。根據(jù)線性代數(shù)理論基本知識有,為了得到頻響函數(shù)和暫態(tài)項的唯一辨識結(jié)果,需要補充2N個復方程。因在任意兩個相異頻率處的頻響函數(shù)之間具有如下關系式:G0(ejωk+l)-G0(ejωk)=∞∑t=0g0te-jωk+lt-∞∑t=0g0te-jωkt=∞∑t=0g0t(e-jωk+lt-e-jωkt)=∞∑t=1g0tφt(ωk+l?ωk)φt(ωk+l?ωk)=e-jωk+lt-e-jωkt}(12)將式(12)代入到式(11)中可得:YN(ωk+l)=G0(ejωk)UN(ωk+l)+VN(ωk+l)+∞∑t=1g0tφt(ωk+l?ωk)UΝ(ωk+l)+1√ΝΝ-1∑t=0ytratejωk+lt(13)式(13)展開可見增加了方程的個數(shù),因?qū)γ恳粋€頻率點k,l的取值為l=-L,…,L。即對每個固定的頻率點k,都存在2L個相鄰點,從而也就使得原來的N個復方程增加到(2L+1)N個復方程。將式(13)中的脈沖響應系數(shù){g0t}∞n+1加入到辨識中,使得有限數(shù)據(jù)中需辨識無窮多個未知參數(shù)。這是無法實現(xiàn)的,為此需降低辨識估計參數(shù)的個數(shù)。假設當時間趨于無窮時,脈沖響應系數(shù)將會衰減至0,即有:g0t=CAtB→0?t→∞ytrat=CAt(x0-xΝ)→0?t→∞}(14)通過截斷(13)式中間兩和項來降低辨識參數(shù)的個數(shù),使得如下的近似式成立:1√ΝΝ-1∑t=0ytratejωk+lt≈1√Νn1-1∑t=0ytratejωk+lt∞∑t=1g0tφt(ωk+l?ωk)≈n2∑t=1g0tφt(ωk+l?ωk)}(15)式(15)中僅采用了n1個暫態(tài)系數(shù),n2個脈沖響應系數(shù),將式(15)代入式(13)中進行截斷可得:?YΝ(ωk+l)=G0(ejωk)UΝ(ωk+l)+1√Νn1-1∑t=0ytrate-jωk+lt+n2∑t=1g0tφt(ωk+l?ωk)UΝ(ωk+l)(16)聯(lián)合n1個暫態(tài)系數(shù)和n2個脈沖響應系數(shù)為整體的未知參數(shù)矢量θ:θ=[ytra0?ytran1-1g01?g0n2]Τ(17)將(16)式改寫成線性回歸模型的形式為:?YΝ(ωk+l)=G0(ejωk)UΝ(ωk+l)+θΤ[1√Ν1√Νe-jωk+l?1√Νe-jωk+l(n1-1)φ1(ωk+l?ωk)UΝ(ωk+l)?φn2(ωk+l?ωk)UΝ(ωk+l)](18)對于式(18)中頻響函數(shù)G(ejωk)和未知參數(shù)矢量θ的求解,可通過如下線性最小二乘優(yōu)化問題來實現(xiàn):argmin{G(ejωk)}Ν-1k=0?θΝ-1∑k=0L∑l=-L|YΝ(ωk+l)-?YΝ(ωk+l)|2(19)對于式(19)的線性最小二乘優(yōu)化問題的求解,可采用經(jīng)典的遞推輔助變量法、最速下降法或高斯牛頓法等。鑒于式(16)的最小二乘問題中有(2L+1)N個復方程,(N+n1+n2)個未知參數(shù)待求量,進一步觀察式(16)中各個量的結(jié)構形式,可采用如下的分離法來實現(xiàn)頻響函數(shù)和未知參數(shù)矢量的求解。4頻響函數(shù)估計的非參數(shù)辨識法對式(16)整理可改寫成如下的特殊矩陣結(jié)構:Y=Φ1G+Φ2θ(20)(20)式中的各個矩陣量的具體結(jié)構為:Y=[Y0?YΝ-1]?Yk=[YΝ(ωk-L)?YΝ(ωk+L)]?G=[G0(ejω0)?G0(ejωΝ-1)]θ=[ytra0?ytran1-1g01?g0n2]?Φ1=[U00?00U1?00??000?UΝ-1]UΝ=[UΝ(ωk-L)?UΝ(ωk+L)]??k?1=[1e-jωk-L?e-jωk-L(n1-1)1e-jωk-L+1?e-jωk-L+1(n1-1)????1e-jωk+L?e-jωk+L(n1-1)]φk?2=[φ1(ωk-L?ωk)UΝ(ωk-L)?φn2(ωk-L?ωk)UΝ(ωk-L)???φ1(ωk+L?ωk)UΝ(ωk+L)?φn2(ωk+L?ωk)UΝ(ωk+L)]Φ2=[?0?1?0?2???Ν-1?1?Ν-1?2]}(21)對于式(20)使用可分離方法的思想在于,第一步是先辨識未知參數(shù)矢量θ,第二步是辨識頻響函數(shù)G。做矩陣Φ2在Φ1所張成空間上的正交投影分解式有:Φ2=Φ‖2+Φ⊥2(22)其中的Φ⊥2分別正交于Φ2和Φ1,根據(jù)矩陣論中的相關性質(zhì)有:Φ⊥2=(I-Φ1(ΦT1Φ1)-1ΦT1)Φ2=diag(I-Φ1(ΦT1Φ1)-1ΦT1)Φ2(23)將式(22)代入到式(20)中可得:Y=Φ1G+Φ‖2θ+Φ⊥2θ(24)在上式的左右兩邊同時乘以(Φ⊥2)*可得:(Φ⊥2)TY=(Φ⊥2)TΦ1G+(Φ⊥2)TΦ‖2θ+(Φ⊥2)TΦ⊥2θ=(Φ⊥2)TΦ⊥2θ(25)由上式可立即得到未知參數(shù)矢量θ的估計值為:?θ=((Φ⊥2)*Φ⊥2)-1(Φ⊥2)*Y(26)將式(26)代回到式(20)中可得:Y=Φ1G+Φ2?θ采用矩陣Φ1的偽逆矩陣可得,頻響函數(shù)G的估計值?G為:?G=(Φ*1Φ1)-1Φ*1(Y-Φ2?θ)=diag{(U*kUk)-1U*k}(Y-Φ2?θ)=[(U*0U0)-1U*0ΔY0?(U*Ν-1UΝ-1)-1U*Ν-1ΔYΝ-1](27)其中:ΔYi,i=0,…,N-1分別為:ΔY0=Y0-[?0?1?0?2]?θ?ΔYΝ-1=YΝ-1-[?Ν-1?1?Ν-1?2]?θ(28)對上述推導過程可歸納為如下頻響函數(shù)估計的非參數(shù)辨識法:(1)采集N對輸入輸出觀測數(shù)據(jù)序列對{u(n),y(n)}Nn=1;(2)根據(jù)式(23)計算由輸入輸出觀測數(shù)據(jù)序列對構成的正交投影矩陣Φ⊥2;(3)根據(jù)式(26)計算增加的未知參數(shù)矢量的估計值?θ;(4)將?θ代入式(28)計算各個ΔYi,i=0,…,N-1;(5)根據(jù)式(27)計算頻響函數(shù)估計值?G。5設置目標函數(shù)的仿真現(xiàn)以某型飛行器的顫振試飛試驗數(shù)據(jù)為例,驗證本文算法的有效性。采用赫伯爾特建立的二維機翼的顫振數(shù)學模型,輸入為人工施加的激勵信號,在整個辨識試驗中,輸入信號可取為常見的偽隨機二元序列信號。輸出是從測點集采集的加速度計測量,時間采樣個數(shù)N=4096。將輸出和輸入的4096個采樣數(shù)據(jù)平均分成4組相等的數(shù)據(jù)塊,每個數(shù)據(jù)塊含有1024個采樣數(shù)據(jù)。真實的模型參數(shù)狀態(tài)空間矩陣分別為:A=[0.65410.6500-0.0009-0.00420.00030.69770.65120.00040.0031-0.0001-0.00070.0050-0.0211-1.01060.01311000001000]B=[-0.00180.001-0.001400]ΤC=[0.1735-0.0301-0.25690.09220.1920](29)從矩陣A的結(jié)構可見,本文顫振試飛試驗的辨識自由度為5。將本文的方法