频域响应辨识优化方法研究

李雅静, 刘旭华, 贾晓鹏, 焦岗

(1.中国飞行试验研究院 飞机所, 陕西 西安 710089;2.西安飞豹科技发展有限公司 设计部, 陕西 西安 710089)

0 引言

等效系统在飞机动态特性评价中是很重要的方法,对试验数据进行时域、频域转换是进行低阶等效系统拟配的基础,试飞数据频域响应分析的质量直接影响拟配效果。

使用快速傅里叶方法进行频域分析时存在两方面问题:第一,需要人工选择窗口宽度[1],且由于单一窗口尺寸不能在整个频率范围内产生满意的相干性,因此单个窗口宽度存在局限性;第二,频域辨识更适用于长时间数据[2],而试飞中使用脉冲、倍脉冲等短时输入较多,包含的频段较窄导致频域辨识结果较差,降低了采集数据的可用性。

本文针对上述问题提出了频域响应辨识的优化方法,即组合窗法和频域辨识修正法。仿真和真实飞行数据的计算结果说明这两种方法能较好地解决前述问题。

1 组合窗法

1.1 窗口选择

首先给出两个公式说明窗口宽度与辨识最小频率和频率响应中的随机误差的关系。

(1)窗口宽度直接决定傅里叶变换中的最小频率fmin:

fmin=1/Twin

(1)

式中,Twin为窗口宽度。

(2)频率响应估计中的随机误差计算公式:

(2)

式中,Cε为常数,用于反应窗口重叠效应,增加窗口重叠比能减小随机误差,代价是快速增加计算量;γxy为相干函数,飞行试验应在风和湍流扰动最小时进行,以获得数据最大信噪比,使相干函数最大化;nd=Trec/Twin为用于求平均的独立时间历程数据个数。较短的窗口Twin能增加nd,但会损害低频范围的辨识精度。式(1)和式(2)表明,大窗口可以增加低频的数据信息、加宽动态范围,而小窗口有助于抑制高频部分的随机误差,所以在不同的频率点应该选择不同的窗口。

1.2 组合窗法

选择单一窗口宽度需要进行折中处理,即为了提高部分频率点的频率响应精度,必须在其他频率点做出牺牲。即使人工找到一个窗口宽度,也仅对所关心的频率范围内有限频段具有最优性。因此需要一种方法能够去除重复的、人工的窗口尺寸优化工作,在所关心的频率范围内生成一个单一的、较准确的频率响应。为此,本文提出了组合窗法。

根据不同窗口宽度在不同频段作用,对动作段采用4个常用窗口宽度32 s,16 s,8 s和4 s进行频域辨识,给出窗口权值曲线如图1所示。32 s窗口有效频段为0.1～5.0 rad/s,16 s窗口有效频段为1～8 rad/s,8 s窗口有效频段为3～10 rad/s,4 s窗口有效频段为5～50 rad/s。

图1 窗口宽度权重Fig.1 Window size weighting coefficient

本文以幅值为例对每个频率点进行加权,相位、相干函数类似,不再陈述。

(3)

式中,GDBC为综合后的幅值;GDBi为相应窗口宽度对应的幅值;γi为相应窗口宽度对应的相干系数;Wi为相应窗口宽度对应的权重值。根据上述加权方法可得到最终的估计幅值、相位和相关函数。以波音747模型为例进行仿真,图2和图3为模型的输入、输出时间历程;根据数据估计得到不同窗口宽度的频域响应如图4所示。

图2 仿真输入时间历程Fig.2 Simulation input time history

图3 仿真输出时间历程Fig.3 Simulation output time history

图4 频域响应Fig.4 Frequency response

图4中r2为相干系数。由仿真结果可以看出,使用小窗口计算得到的频率响应在低频段相干性较差,使用大窗口计算得到的结果在高频段随机误差较大,而使用组合窗法计算得到的结果在整个频段都有较好的相干性且较为平滑。辨识结果与建模传递函数频响比较如图5所示。可以看出二者十分接近,表明使用组合窗法不需人工选择窗口宽度,并能得到一个单一的、在整个频率范围内较精确的频域响应。

图5 真实频响与辨识结果比较Fig.5 Comparison between true frequency response and identification result

2 频域辨识修正法

频域辨识方法适合长时间数据,时域辨识方法适用于短时间数据。试飞中脉冲、倍脉冲等短时间动作较多、数据较短,对时域辨识结果进行频域分析,既可用于修正频域辨识结果,也可用于验证频域辨识结果,最终能得到较好的辨识结果。

2.1 时域辨识方法

设飞机的高阶增稳系统及其等效系统在相同的输入u(t)作用下的输出响应分别为yH(t)和yL(t),则两个系统近似为输入、输出等效的条件。这里yL(t)并不是一般意义上的低阶等效系统,而是一个若干阶的系统。

代价函数[3]为:

(4)

yH(t)可以通过直接测量得到;yL(t)在第k个采样时刻的数值应满足如下差分方程:

A(q-1)yL(k)=B(q-1)u(k)+e(k)

(5)

其中:

A(q-1)=1+a1q-1+…+anq-n

(6)

B(q-1)=b1q-1+…+bnq-n+1

(7)

式中,e(k)为方程残差,由剩余高阶特性、随机风干扰及测量噪声等引起。

MATLAB自带的系统辨识工具箱[4]使得辨识方法十分简便,上述时域辨识方法只是其中的一种ARX,还包括PEM,BJ,OE,ARMAX等。时域法辨识采用多个模型结构和不同的辨识算法。使用不同辨识函数进行辨识后,比较辨识系统输出与真实系统输出的相似度并选出最优辨识结果,即选择相似度最大的一组进行频域响应计算,得到时域辨识系统的频响。

2.2 时域辨识对频域辨识修正法

使用时域方法得到的频响修正直接频域辨识结果。修正方法为:分别对数据进行时域辨识和频域辨识,时域辨识得到状态空间方程或传递函数,再对该传递函数进行频域分析,得到系统频域响应;频域辨识直接使用快速傅里叶变换进行频域分析;将两个结果进行比较并加权得到最终的频响。

首先使用杆输入激励由系统辨识工具箱得到的时域模型,若输出与飞机真实输出相似度大于85%,则说明该模型可以代替飞机真实模型;若相似度小于85%则不修正。修正量计算式为:

(8)

(9)

Wf(i)=1-Wt(i)

(10)

修正公式为:

G(i)=Wt(i)Gt(i)+Wf(i)Gf(i)

(11)

P(i)=Wt(i)Pt(i)+Wf(i)Pf(i)

(12)

式中,r2(i)为频域辨识的相干系数;fitmax为时域辨识的相似度;Wt(i),Wf(i)分别为时域辨识频响与频域辨识频响所占比例;Gt(i),Gf(i)分别为各频率点时域辨识与频域辨识幅值响应;Pt(i),Pf(i)分别为各频率点时域辨识与频域辨识相位响应。由式(9)可见，加权系数由频域辨识的相干系数和时域辨识的相似度组成,哪个大,则响应辨识方法的结果所占权重大。

一般情况下,使用相干函数的大小评价频域计算结果,图6为某机俯仰角速度频域辨识结果。直接辨识结果的低频段相干函数小,拟配时一般适当缩小拟配频率范围,缩小到0.7～10.0 rad/s[5]。图7所示的拟配结果并不好,因为该机的频率范围较低,从0.7 rad/s开始拟配频率范围较窄,导致拟配结果差。

图6 俯仰角速度频域辨识结果Fig.6 Frequency domain identification results of rate of pitch

图7 时域修正前俯仰角速度拟配结果Fig.7 Matching results of rate of pitch before time domain

在低频段辨识结果差、缩小拟配范围也无效的情况下,使用时域结果或综合结果可以得到较好的拟配结果,如图8所示。大量试飞数据计算结果表明,时域辨识对频域辨识的修正方法大大提高了试飞数据的利用率。此外时域结果可以对频域结果进行验证,如图9所示。在计算过程中,如果时域估计模型的频响与直接计算得到的频响较接近,一般情况下可以得到很好的拟配结果。

图8 时域修正后俯仰角速度拟配结果Fig.8 Matching results of rate of pitch after time domain

图9 俯仰角速度频响Fig.9 Frequency response identification result of rate of pitch

3 结束语

本文从实际出发,提出了组合窗法、时域辨识对频域辨识修正的方法,很好地解决了频域响应辨识过程中存在的问题,通过实例验证了方法的有效性。将本文方法应用于实际飞机动态特性指标计算,得到了很好的效果,使辨识结果更准确、拟配结果更合理,而且大大提高了试飞数据的利用率。

参考文献：

[1] Tichler M B,Remple R K.Aircraft and rotorcraft system identification[M].USA:American Institute of Aeronautics and Astronautics,2006:145-167.

[2] Klyde D H,Bachelder E N,Thompson P M,et al.Flying qualities parameter identification using short duration flight test inputs[R].AIAA-2007-6385,2007.

[3] 马维金,王俊元,李凤兰,等.基于ARX模型的控制系统辨识及稳定性分析[J].中北大学学报,2010,31(1):9-13.

[4] 徐昕,李涛,伯晓晨.MATLAB工具箱应用指南:控制工程篇[M].北京:电子工业出版社,2000:28-41.

[5] 包立平,李春锦.评定增稳飞机飞行品质的纵向频域等效系统法[J].北京航空学院学报,1984,(4):31-43.