基于离散傅里叶变换的频响函数非参数辨识

基于离散傅里叶变换的频响函数非参数辨识

频响函数是振动结构中输出信号和输入信号的傅里叶变换比率。振动结构的频响函数包括原振动结构物理参数的所有信息,如模态参数的确定、振动结构损失的程度。反之振动结构物理参数的信息也唯一地决定了频响函数的变化。目前的研究侧重直接从频响函数中辨识振动结构的模态参数,相对于模态参数而言,频响函数更直接、误差小、含有更丰富的原始数据信息。正是频响函数的这些优势,可将频响函数作为模态参数辨识和振动结构监测损伤的主要依据。虽对频响函数的研究较广泛,但绝大部分的研究都是假设已知的频响函数事先存在,或利用原始定义直接得到频响函数。文献为了克服对输入-输出观测序列作离散傅里叶变化时所带来的暂态和泄露谱影响,将通常的离散傅里叶变化推广到同时考虑初始和终端状态。采用频响函数数据的频域多输入多输出状态空间模型,用频响函数数据代替传统的状态空间模型中的输出,进而采用子空间辨识法求解状态空间模型中的系统矩阵。但文献中的子空间辨识法能使用的前提条件是在各个频率点处的频响函数估计值在统计概率意义下达到渐近有效性和一致性。文献分析振动结构在环境激励下的模态参数辨识问题,其使用的辨识方法仍是子空间辨识法,以绕开频响函数的不准确而使得直接由不准确的频响函数得到较大偏差的模态参数值。文献直接利用频响函数识别结构损伤,即频响函数曲率法的原理类似于振型曲率法,但不需要测试振型,比振型曲率法识别效果更佳。以上代表在振动工程领域对频响函数的研究还仅仅停留在频响函数的应用价值取向上,还未对频响函数的本质做进一步的理论研究。但在系统辨识领域已对频响函数做了一些基础性的铺垫工作。如文献提出了多种不用于原始定义的频响函数估计表达式,并分析了这些表达式各自的统计性质。文献从二阶统计量-方差的角度分析为得到渐近无偏的频响函数估计值,需要在对输入-输出观测序列的离散傅里叶变换过程中引入加权滤波器,以达到平滑滤波噪声的功能。文献分析线性动力系统中频响函数的频域极大似然估计。文献分析非线性系统的最优线性无偏估计,即用一个频响函数估计值来最优地逼近于原非线性系统。本文深入研究频响函数的辨识估计,为振动工程领域中频响函数的直接使用提供准确表达式,将频响函数估计方法应用于振动结构的模态参数辨识之中。频响函数的估计通常有两种方法:参数法和非参数法。频响函数估计的参数法中,当增加了额外的未知参数势必造成参数辨识成为一个奇异问题,不能得到关于未知参数的唯一解。因此本文研究频响函数估计的非参数法。根据频响函数的定义可知,首要一步是对输入-输出观测序列做离散傅里叶变换,在此变换过程中同样考虑初始和终端状态带来的暂态泄露项和观测噪声谱项对频响函数估计的影响。在增加此两项的基础上,为得到频响函数估计值,联合频响函数、初始-终端状态和脉冲响应系数作为整体的待辨识未知参数矢量,将频响函数估计问题转化为求解一个线性最小二乘优化问题。针对此线性最小二乘优化问题的特殊形式,提出一种可分离的求解过程1结构模态参数辨识频响函数主要应用于模态参数辨识中,在时域环境下进行模态参数辨识,通常采用微分或差分方程、传递函数矩阵等形式来描述系统结构的振动行为。而以上输入-输出模型只刻划了结构的外部特性,并未深入内部。而状态空间模型则可以深入到结构的内部情况。此处介绍如何将输入输出模型转化为状态空间方程模型,以及得到的频响函数与状态方程模型之间的关联。在振动工程领域经常遇到振动响应问题:对于非时变、粘性阻尼、线性自由度为N的系统,可以根据Hamilton原理,推导出在外力u(t)作用下的运动方程为:Μ⋅⋅y(t)+C1˙y(t)+Κy(t)=u(t)(1)式中:M,C,K∈RN×N分别为振动系统的质量矩阵、粘性矩阵和刚度矩阵,y(t)是位移向量。式(1)给出系统的广义运动方程,且该式是结构动力学研究最常用的形式,也是结构模态参数辨识的基本方程式。将上式与某些恒等式组合在一起可得:˙y(t)=Ι˙y(t)⋅⋅y(t)=-Μ-1Κy(t)-Μ-1C1˙y(t)+Μ-1u(t)}(2)假设状态变量x(t)为2N维矢量:x(t)=[y(t)˙y(t)]Τ可得状态方程为:˙x(t)=[˙y(t)⋅⋅y(t)]=A1x(t)+B1u(t)(3)其中:A1,B1分别为:A1=[0Ι-Μ-1Κ-Μ-1C1]‚B1=0Μ-1Ι其中:A1∈R2N×2N为系统矩阵,B1∈R2N×1为输入矩阵。系统的输出向量与状态向量x(t)间的关系为:y(t)=Cx(t)+v(t)(4)其中:y(t)为m维的输出向量,C为m×2N维的结构输出矩阵,v(t)为m维的系统测试输出噪声。联合(3)、(4)式所述状态方程的离散化形式为:xn+1=Axn+Bunyn=Cxn+vn}(5)通常对(5)式进行离散傅里叶变换可得,此时经常忽略初始状态和终端状态下,系统的频响函数矩阵为:G(z)=C(zI-A)-1B(6)求取状态矩阵A的特征值和特征向量,则有如下等式成立:A=PΛP式中:Ρ=[φ(1)φ(2)⋯φ(2Ν)]Τ,且φ(i)为A的特征向量,Λ=diag(λ1λ2⋯λΝ),且λi为A的特征值。代入可将状态空间模型转化为模态展开模型:G(z)=CΡ(zΙ-Λ)-1Ρ-1B=[φφ*](zΙ-Λ)-1[LL*](7)上式中称矩阵A的特征值为系统的极点,CP为模态的振型矩阵,φ-1P为模态的贡献因子矩阵。以上的具体分析在于系统矩阵A的子空间辨识求解,但若能求得频响函数G(z)的估计值ˆG(z),则可以利用频响函数估计值ˆG(z)来辨识模态参数。频响函数矩阵式(6)中第(i,j)个元素可表示在第j自由度方向施加单位简谐荷载时,在第i自由度方向产生的位移响应。响应测点i和激励点j之间的位移频响函数为:Gij=Ν∑r=1φirφjrΚr-ω2rΜr+jωrCr(8)式中:ωr为振动结构的第r阶模态频率;Kr,Mr,Cr是对应的模态刚度、模态质量和模态阻尼;φir,φjr分别是对应的(i,j)点处的振型幅值。由式(8)可知,若左边的频响函数估计值求出后,可将其改写成分式和的形式,从各个分式中的系数即可直接得到对应各个模态参数值。2散傅里叶变换重写式(5)所示的离散化状态空间模型形式:xn+1=Axn+Bun,yn=Cxn+vn上式的频响函数定义为:G0(ejω)=∞∑t=0g0te-jωt=∞∑t=1CAtBe-jωt(9)其中:g0t=CAtB表示系统的脉冲响应系数,本文的研究目标在于,从一族输入-输出观测序列对{u(n),y(n)}Νn=1(N表示观测数据个数)中,辨识出在所有频率网格点ωk=2πΝk‚(k=0⋯Ν-1)处的频响函数G0(ejωk)的非参数估计值ˆG0(ejωk)。当N个输入-输出观测数据对可获得时,应用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的离散傅里叶变换过程中,考虑系统的初始状态x0和终端状态xN的存在,且推广了式(6)中的一般广义情况。将式(10)应用于式(5),并经过整理可得: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为初始状态,xN为终端状态,第2项TN(ejωk)表示暂态项或泄漏谱,其产生的原因在于有限观测时间下初始和终端状态带来的影响。且由TN(ejωk)的表达式可见,只有当用周期信号激励系统时,有x0=xN,从而暂态项可消除。对于式(11)中的真实频响函数G0(ejωk)的估计值求法通常采用如下定义等式:ˆG(ejωk)=YΝ(ωk)UΝ(ωk)+ΤΝ(ejωk)UΝ(ωk)+VΝ(ωk)UΝ(ωk)此估计值的特点在于,估计值ˆG(ejωk)仅定义在固定的频率点处;此估计值是无偏估计,且其方差随着1/N而递减;因没有平滑滤波的功能存在,该估计值的方差不会随着N的增加而递减;在随机激励作用下,该估计误差将会任意大。因此根据输出信号和输入信号的离散傅里叶变化之比所得到的频响函数估计是不理想的,不能直接应用于振动模态参数辨识之中。为了从输入输出观测数据序列的离散傅里叶变换UN(ωk)、YN(ωk)中辨识在整个频域范围ωk=2πΝk‚(k=0⋯Ν-1)中的频响函数G0(ejωk)的估计值ˆG(ejωk)。3线性关系的建立因式(11)中不仅含有频响函数,还含有暂态项,同时将暂态项和频响函数作为需要辨识的量。但式(11)展开后仅有N个复方程(k=0,…,N-1),而辨识过程中却有3N个未知待辨识量。根据线性代数理论基本知识有,为了得到频响函数和暂态项的唯一辨识结果,需要补充2N个复方程。因在任意两个相异频率处的频响函数之间具有如下关系式: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)展开可见增加了方程的个数,因对每一个频率点k,l的取值为l=-L,…,L。即对每个固定的频率点k,都存在2L个相邻点,从而也就使得原来的N个复方程增加到(2L+1)N个复方程。将式(13)中的脉冲响应系数{g0t}∞n+1加入到辨识中,使得有限数据中需辨识无穷多个未知参数。这是无法实现的,为此需降低辨识估计参数的个数。假设当时间趋于无穷时,脉冲响应系数将会衰减至0,即有:g0t=CAtB→0‚t→∞ytrat=CAt(x0-xΝ)→0‚t→∞}(14)通过截断(13)式中间两和项来降低辨识参数的个数,使得如下的近似式成立: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个暂态系数,n2个脉冲响应系数,将式(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)联合n1个暂态系数和n2个脉冲响应系数为整体的未知参数矢量θ:θ=[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)中频响函数G(ejωk)和未知参数矢量θ的求解,可通过如下线性最小二乘优化问题来实现:argmin{G(ejωk)}Ν-1k=0‚θΝ-1∑k=0L∑l=-L|YΝ(ωk+l)-ˆYΝ(ωk+l)|2(19)对于式(19)的线性最小二乘优化问题的求解,可采用经典的递推辅助变量法、最速下降法或高斯牛顿法等。鉴于式(16)的最小二乘问题中有(2L+1)N个复方程,(N+n1+n2)个未知参数待求量,进一步观察式(16)中各个量的结构形式,可采用如下的分离法来实现频响函数和未知参数矢量的求解。4频响函数估计的非参数辨识法对式(16)整理可改写成如下的特殊矩阵结构:Y=Φ1G+Φ2θ(20)(20)式中的各个矩阵量的具体结构为: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)使用可分离方法的思想在于,第一步是先辨识未知参数矢量θ,第二步是辨识频响函数G。做矩阵Φ2在Φ1所张成空间上的正交投影分解式有:Φ2=Φ‖2+Φ⊥2(22)其中的Φ⊥2分别正交于Φ2和Φ1,根据矩阵论中的相关性质有:Φ⊥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)由上式可立即得到未知参数矢量θ的估计值为:ˆθ=((Φ⊥2)*Φ⊥2)-1(Φ⊥2)*Y(26)将式(26)代回到式(20)中可得:Y=Φ1G+Φ2ˆθ采用矩阵Φ1的伪逆矩阵可得,频响函数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)对上述推导过程可归纳为如下频响函数估计的非参数辨识法:(1)采集N对输入输出观测数据序列对{u(n),y(n)}Nn=1;(2)根据式(23)计算由输入输出观测数据序列对构成的正交投影矩阵Φ⊥2;(3)根据式(26)计算增加的未知参数矢量的估计值ˆθ;(4)将ˆθ代入式(28)计算各个ΔYi,i=0,…,N-1;(5)根据式(27)计算频响函数估计值ˆG。5设置目标函数的仿真现以某型飞行器的颤振试飞试验数据为例,验证本文算法的有效性。采用赫伯尔特建立的二维机翼的颤振数学模型,输入为人工施加的激励信号,在整个辨识试验中,输入信号可取为常见的伪随机二元序列信号。输出是从测点集采集的加速度计测量,时间采样个数N=4096。将输出和输入的4096个采样数据平均分成4组相等的数据块,每个数据块含有1024个采样数据。真实的模型参数状态空间矩阵分别为: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的结构可见,本文颤振试飞试验的辨识自由度为5。将本文的方法